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ABSTRACT 


A method,  made  up  of  several  algorithms,  is  given  for  computing  the  incomplete 
gamma  function  ratio,  P(a,  x)  and  its  complement  Q(a,  x)  for  all  real  arguments 
a > 0,  x 5*  0.  The  difficult  case  when  a and  x are  large  is  treated  by  a modified  version 
of  Takenaga’s  method.  The  resulting  computer  program  is  efficient,  yields  both  P and 
Q correctly  to  within  1 unit  in  the  twelfth  significant  digit,  or,  at  the  user’s  option,  to 
within  one  unit  in  the  sixth  or  third  significant  digit.  A FORTRAN  listing  of  the  pro- 
gram is  included. 
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1 INTRODUCTION 


The  object  of  this  report  is  to  describe  an  efficient  procedure  for  computing  the 
incomplete  gamma  function  ratio 


P(a,  x)  = >(a,  x)/P(a).  a > 0,  x > 0, 

and  its  complement 

Q(a,  x)  = l'(a,  x)/l'(a),  a > 0.  x>0, 

to  within  a unit  in  twelve,  six,  or  three  significant  digits,  where 


7 ( a , x ) 


■X** 


ta  1 dt. 


(1) 


f2) 


(3) 


l 


' ( a , x ) = / e 1 t J 1 d t . 


><a,~)  = I (a,  t))=  1(a)  = / e ’ t dt 

'0 


-f 


(4) 


(5) 


with 


0 **  l‘  • 1 . 0 • 0 = I i»  v;  | . 

The  computing  program,  in  I OUTRAN.  designed  for  the  CIX-6700*.  uses  7 
different  mathematical  formulations  for  I’  ami  ().  the  asymptotic  expansion  (a  • °°)  for 
the  complete  gamma  function.  1(a).  and  a series  expansion  for  its  reciprocal.  They  are 
all  given  in  the  next  section  with  their  domains  of  application.  The  algorithm  used 
with  each  formulation  is  given  in  Section  5 

If  the  user  desires  1 2 significant  digits  lor  I’.  ().  or  hot U , a parameter  r is  set.  in  the 
call  to  the  program,  to  the  value 

.=(,=5x10  IJ  = 5(  13)  (6) 


*The  ( IX  -b^UO  is  a large  scale  binary  computer  capable  of  one  million  arithmetic 
operations  per  second.  It  has  a 60  bit  binary  word  length  of  which  48  are  used  to 
express  the  mantissa  of  a number. 


I 


— 


If,  on  the  other  hand,  the  user  wants  only  6 or  3 significant  figures,  then  e is  set, 
respectively,  to 


e = e2  = 5 x 10-7  = 5(-7).  e = fj  = 5 x lO'4.  (7) 

Besides  three  available  levels  of  accuracy  for  the  output  quantities  P and  Q,  the 
program  also  has  the  capability  to  determine,  in  most  cases,  if  the  input  quantities  a 
and  x are  contained  in  “zones  of  non-computation.”  These  are  domains  in  the 
ax-plane  for  which  P < e or  Q < e;  tor  sucli  (a,  x).  1*  is  set  to  zero  and  Q is  set  to  one 
or  vice  versa  depending  on  the  values  of  (a,  x).  Section  3 deals  with  the  ways  in  which 
the  program  senses  on  the  above  inequalities. 

The  average  computing  time  per  case,  when  r = t , (b  significant  digit  accuracy  in 
P and  Q)  is  about  1 millisecond:  for  e = <. , (12  significant  digits  in  P and  Q)  the  average 
computer  tirm  per  case  is  about  15  to  20' - longer  and  for  c = c3  about  10  to  15% 
shorter. 

The  validation  of  the  program  is  described  in  Section  5.  A complete  FORTRAN 
IV  listing  is  given  in  Appendix  I . 

There  are  many  applications  for  the  incomplete  gamma  function  ratios  such  as 
computing  the  Poisson  distribution.  | I . p.  ‘>5l>| . Perhaps  their  most  well-known  use  is 
for  computing  the  Chi-square  distributions  P(\:|r)  arid  Q(\*|iO  = I - P(\2  nO.  where 

P(\ 3 |r)  = [2‘  : I'(|'|2)1  1 f e " : u"  1 du.  (8) 

-z(> 

By  letting  u - 2t  in  (8),  the  relationship  between  (\:.  v)  and  (x.  a)  becomes  evident. 


J 0 


e 1 t,,/;)  1 


dt/l  (i  2 ) - PO/2,  \ * 2). 


•More  precisely,  one  should  state  that  in  the  case  of  ((>).  P and  Q are  given  correctly  to 
within  one  unit  in  the  twelfth  significant  digit,  and  to  within  one  unit  in  the  sixth  or 
third  significant  digit  in  the  case  of  (7). 


so  that 


a = v/2,  x = \2/2  . 

The  extensive  use  of  the  Chi-square  distributions  in  statistics  is  sufficient  to  warrant 
a good  machine  program  for  their  equivalent,  the  incomplete  gamma  function  ratios. 
Surprisingly,  to  the  best  of  our  knowledge,  there  is  only  one  computer  program  docu- 
mented in  the  literature,  [11].  In  any  case,  the  authors  are  not  aware  of  a program  for 
P and  Q that  has  the  overall  speed,  versatility,  accuracy,  and  range  of  the  present 
program. 


1 

I 

I 

f 
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2.  MATHEMATICAL  FORMULATIONS  AND  DOMAINS  OF  APPLICATION 


In  this  section  we  summarize  the  mathematical  formulations  and  state  the  domains 
on  which  they  are  used. 

It  is  worth  noting,  since  we  require  relative  accu'acy,  that  it  is  necessary  to  com- 
pute Q(P)  if  P(Q)  is  near  one  (and  then  find  P(Q)  from  1 Q(P)).  Consequently,  the 

program  is  designed  to  compute  : if  i > 1 and  a > x.  because  in  this  case  P cannot  be 
near  one.  Similarly,  if  1 < a < x.  a > in  10.  then  0 is  computed  for  it  cannot  be  near 
one.  If.  however,  a < 1 and  x "4  a then  I*  may  be  r ;ar  one  and  Q should  be  computed, 
although  if  x is  sufficiently  small  P will,  of  course,  be  near  zero  and  should  be  com- 
puted. This  behavior  of  P near  the  origin  of  the  ax-plane  results  from  the  nonexistence 

of  lim  P(a.  x).  a fact  easily  concluded  from  (10)  below, 
x • 0 
a - 0 


I igure  1 shows  the  domains  over  which  the  various  formulations  are  used.  It  and 
the  flowchart  on  page  44  implement  the  remarks  mat  follow.  It  is  assumed  P > t 
and  Q > r . 

We  proceed  to  specify  the  formulations  and  their  domains  of  application. 


K(a.  \ 1 
P - 

' \k 

. M.  p-2631  (10) 

a 

, , (a  + 1 1 (a  + k) 

Rta.  x)  = e ' v’  Tia) 

(ID 

a 1 00 

1 ) 1 a. 

a • x . 

2)  2a  is  not  an  integer  when  \ • a 

a) 

s*  a • \ ,n  10 

h)  a 

<1,  x < 1.5,  !'  s;  0.40 

c)  a 

<1,  x 1.5.  R > 0.101  x ( ( x + 2 a )/( \ + 1)| 

a>  100 

3a  > 4x 

4 


m ;p 

III 

• •a 

* ' 1 

. 

1 

E 

t: 

y — \ 

►4  •-» 

i 

1 1 * — 

K(a.  x)  * (a  ! ) • • (a  k) 

o«».  — - — > + ■£  — 


a<  100 


I < a < x.  x ^ 31(e  = e,  I.  x > 17(e  = e2 ).  x 5*  9.7(c  = )* 


a > 1 00  4x  > 5a 


Q(k  +g.  x)  = Q(k  +g  1 . x)  + R(k  + g 1 , x)/(k  + g 1).  ( 1 . p.  262 J (13) 

j I a is  an  integer 
g=  \ 

M/2  a is  not  an  integer:  k = 1.  2.  • • . a g 

/°°  : 

e u du  = erfc(\/x).  g = 1/2 

Q(g,  x)  = : v'x-  ,I4\ 


g=  I. 


a < 100 


2a  is  an  integer,  a < x < 3 i . 1 7.  9.7  for  e = e( , c, , e3,  respectively. 


Q(a,  x)  = R(a,  x) 


1 10,  p.  356] 


a < 100 


1 . 2a  is  not  an  integer 

a)  1 < a<  x,  fn  10<  x<  31,  17,9.7fore  = €,,€,,€3, 

respectively. 

b)  a<l,x>1.5,R<O.I01x(x  + 2 a)/(x  + 1) 


Q(a,  x)  = J (H  + L)  + J(H  + L)  HL  + JHL 


(16) 


J- 


( x)k 


k _ i (a  + k)  k! 

A (a  fn  x)k 

^ k ' 
k 1 


Z f. 


,k  - I 


H = 


k = 2 


(1/a) 


(1  a)+  Z Ck(a  !)k  ' 


0 < a < 1/2 


1/2  < a < 1 


The  coefficients  C'k  are  given  in  Appendix  A,  also  see  (25).  Discussion  of  ( 1 6) 
given  below  and  on  pages  10  and  1 1 . 

a < 100 

a < 1,  x < 1.5,  P > 0.90. 


The  need  for  (16)  arises  when  a < 1,  since  then  it  is  not  true,  even  though  x < a, 
that  P cannot  be  near  one.  In  fact,  for  a = x s 0.038,  P s 0.90,  for  a = 10~3, 
x = 10-\  P = 0 9913  • • • , and  for  a = 10"5,  x = 10“  P = 0.9997  • • • . Thus,  in 

such  cases  it  is  Q,  rather  than  P,  which  must  be  computed.  This  behavior  of  P near  the 
origin  of  the  ax-plane,  as  mentioned  previously,  can  be  accounted  for  by  noting  from 
(10)  that  the  double  limit  for  P as  x ♦ 0,  a * 0 does  not  exist. 


If  a ^ 100  and  3a  < 4x  < 5a,  none  of  the  above  formulations  suffice  for  comput- 
ing purposes.  (If  4x<  3a,  (To)  is  used  and  if  4x  > 5a  then  (jj)  is  UJed.)  In  this 
main  area  of  difficulty,  we  have  found  that  Takenaga's  method  (9) , with  appropriate 
changes,  works  quite  well.  The  final  formulas  are  given  by  (17#  (24).  The  mathe- 

ma*ical  formalism  is  easily  programmed  into  an  accurate  and  relatively  efficient  proce- 
dure for  computing  P(a.  x)  when  x < a 1/3.  or  Q(a,  x)  when  x > a 1/3.  The 
details  for  deriving  (17)  (24)  are  given  in  Section  4. 


Let 


■ "'TV'  T1V 


r 


T(a,  x)  = 


L(a,  x)  if  x < a - 1/3 

I Q(a,  x)  if  x > a - 1/3. 

Then  Takenaga’s  analysis,  after  some  changes,  reduces  to 


s C 

Bn 

1 

B, 

~~  (B3 

— B ) 

0 

12a2  2 

18a4  l 3 

16  7 

( 1 7)" 


1 / 37  1 \ 

B, B + — BJ 

24a6  \ 4 225  5 432  7 


30uR  'Bs 


743  49 

B.  + B, 


302^'  6 4320  7 82944 


■) 


— (—  B,, 2.3414  67226  51062  x 10‘22  B,J 

a24  '78  ” 24/ 

+ M + _L  U J1A  +±A) 

15a3  3uf  ' 7 60'  27a7  \ 4 140  s 160  7 

I / I 193  1187  i \ 

a'1  '33  22680  2268000  J55520  7 


1 / 1 

— - A„  • • • 2. 

«25  '81  1 


2478  08537  45020  x 10 


*The  numerical  coefficients  in  (17)  and  (18)  were  computed  by  Alfred  H.  Morris 
using  his  "I  lap”  algebraic  routine.  ( 5 i . The  complete  computed  sets  of  coefficients 
are  given  in  Appendix  B. 


where,  with  b = a 1 


1 


31  , 3413  , 361733 

C = 1 + — b- 1 b 2 + b-3+ b-4 

36  2592  1399680  201553920 


(18) 


13888281  . 7565202533 

b 5 + b-6 


50791587840  7836416409600 

a = 3(a  - l/3)l/2 


(19) 


The  Ak  and  Bk,  for  a given  integer  k,  have  different  meanings  depending  on  whether 
x < a - 1/3  or  x > a - 1/3.  If  x < a - 1,3,  then 


1 /~S  2 

Ak  - — — / z2kt|  e~z  12  dz 
\/Tn  J — oo 


(20) 


b''jt  LS  dz' 


k = 0,  1, 


(21) 


where 


s = a 


\a  - 1/3/ 


^ 0. 


(22) 


If  x > a - 1/3,  then  (20)  and  (21 ) arc  replaced  by 

1 


A = — r f z2k*'  e 1 I2  dz 

\Z2ir  J s 


(23) 


I C 00  2 

B = / z2k  e_z  !2  dz, 


k = 0,  1,2, 


(24) 


In  this  case  also,  s is  given  by  (22),  but  now  it  is  positive.  We  remark  that  the  program 
uses  a slightly  different  form  for  (17).  See  (115), 


The  evaluation  of  R(a,  x)  (see  (11)),  which  contains  l/r(a),  is  carried  out  in  two 
different  ways.  If  a < 30,  then  the  factor  x*  e“*  is  computed  in  the  form  e-**t(nx, 
and  subsequently  r(a)  is  evaluated  as  a separate  factor  by  the  recurrence  relation 


r(j  + x)  = (j  + x - l)  r(j  + x - l),  j = l,  2,  • • • , a - x,  x # o. 


. -iriVriirril  -n 


•'  ’ 'll  !!.■  , ,,  > i,  F.„, 


" '"I’l.1  (11,1^1.  1.1  , J1  I 


where  X denotes  the  fractional  part  of  a and  a > 1.  The  procedure  is  initiated  by  com- 
puting l'(X)  from  a series  expansion  for  its  reciprocal.  The  quantity  1/P(X)  is  computed 
to  14  significant  digits  by  the  polynomial  approximations,  [131, 


17 


X Itjc,  X1  ' 


k 2 


l/r<x)  = 


l + £ck(x  l)k  ' 


k 2 


0<  X < 1/2 


1 12  < X < 1 ; 


(25) 


the  Ck's  are  given  in  Appendix  A.  By  using  the  first  equation  of  (25)  and  P(X)  = 
T(1  + y)  = y P(y)  where  X = 1 + y.  so  that  1/2  < X < 1 , 1/2  < y < 0,  the  second  equa- 
tion of  (25'  is  obtained  If  X = 0,  then  the  above  recurrence  relation  is  started  with 
j = 2 (instead  of  one)  and  l’(  1 ) = I. 

If  30  < a,  then  the  asymptotic  expansion  for  fn  P(a),  (a  ♦ -).  [1,  P-  257), 


i'nP(a)2(a  1,2)  fn  a a + - fn  2n  + J.(a), 


(26) 


where 


I 

L(  a ) = — + 


1 


I 


- • 0<  fl  < ij  (27) 


12a  360a3  1260a5  1680a1 

is  combined  with  the  logarithm  of  the  numerator  of  ( 1 1 ) to  give,  after  exponentiation. 


I 


R(a,  x)  s CXp 

s/2ff 


(a  x)  + a fn L(a) 

a 


(28) 


It  is  easy  to  show  that  the  argument  of  the  exponential  in  (28)  is  always  negative  for 
admissible  values  of  a and  x.  Hence  no  scaling  problems  can  occur  due  to  the  large 
magnitude  of  a and  x. 

We  do  not  give  detailed  derivations  of  the  above  equations  except  for  (17)-  (24) 
which  are  developed  in  Section  4.  Equations  (10),  (12),  and  (13)  are  easily  obtained 
using  integration  by  parts  on  (1 ) or  (2);  (14)  follows  directly  from  (2)  when  a = 1/2  and 
t = uJ ; ( 1 5)  is  derived  in  [10].  Equation  (16)  follows  by  substituting 


10 


♦ . 

t.i  nh  «!«■■  


I 1}  IIWPL  ^9>iU  ' '■ 


1 


p = 

xJ/l'(a  + 1 ) 

1+af  ("X)k 

k = 1 (a  + k)  k! 

(See  ( 1 , p.  262) ) 


x*  = 1 + £ 

k - I 


(a  fn  x)k 
k! 


(29) 


(30) 


1 /r  (a  + 1)  = 


0<  a<  1/2,  (See  (25)) 


(31) 


1/2  < a < I. 


in  the  expression  1 - P(=Q)  The  derivation  of  the  series  expansion  for  1/P(a)  is 
sketched  in  1 1 3 1 . The  asymptotic  expansion  lor  i n l‘(a)  is  discussed  in  |6,  p.  2931. 


Table  1 has  been  compiled  to  give  some  ide.i  of  the  number  of  terms  of  (10), 
(12).  or  (15)  that  are  needed  for  various  arguments  (a.  x)  with  t = r , and  c = e2 . The 
third  column  of  the  table  identifies  the  equation  that  was  used  by  its  r,'mber  in  the 
text.  The  fourth  and  fifth  columns  list  the  number  of  terms  of  that  equation  needed 
for  e = and  c = e , respectively.  For  example  a = '/.I . x = 28  uses  1 2 terms  of  ( 1 5) 
forc(  and  7 terms  of  ( 12)  for  c . No  examples  are  given  for ; 1 3)  since  the  number  of 
iterations  is  always  the  integer  part  of  a.  Also  no  data  is  given  for  ( 16).  Since  ( 16)  is 
used  when  a < 1 and  x < 1 .5,  it  yields  an  efficient  algorithm.  In  the  worse  case,  when 
x - 1.5  and  r = c ( . no  more  than  18  terms  of  J *re  needed:  the  evaluation  of  L requires 
less  then  12  terms  and  in  the  worse  case  of  a 2 1/2,  no  more  than  17  terms  of  II  are 
used. 


In  the  next  section  we  show  how,  for  small  c > 0.  advantage  is  taken  of  a property 
of  I*  and  Q.  Namely,  if  a is  large,  then  x cannot  be  too  far  from  a.  otherwise  P < c or 
I P = Q < e depending  on  whether  a > x,  or  x > a.  respectively.  In  other  words, 
there  is  a very  limited  region  of  the  • x-planc  for  whi<.\  c < P < I c and  for  which 
extensive  calculations,  using  one  of  the  above  formulations,  is  necessary.  Inequalities 
will  be  derived  which,  in  the  case  of  equality,  give  very  close  approximations  to  the 
boundaries  of  this  region.  Farlier.  we  referred  to  the  exterior  and  boundary  of  this 
region  as  the  “zone  of  non-computation."  The  computing  time  for  any  one  of  the 
inequalities  is  relatively  short.  As  a result,  the  average  computing  time  per  case  is  mark- 
edly reduced,  since  the  average  computation  time  for  computing  P and  Q.  when  (a,  x) 
is  in  the  “zone  of  computation”  (c  < P < I t ) is  generally  much  longer  than  when  it 
is  not. 


| 


i 


I I 


j 


V 


5 


10  I 

i|  15 
- ^ 

I 1 

2.1  I 

7 1 

13  I 

]TjJ 

25  JJ. 

32  i; 


20  10 

30  10 

32  12 

50  12 


29  I ‘l 

16  10 


99.99  50 

75 
99.9 
100 


20  I 


28  I 1 


185 


1 2/^ 

/ * 

20 

* 

10/ 

/ * 

35 

* 

10 

56 

35 

10 

86 

61 

12 

65 

49 

12 

54 

37 

12 

43 

27 

i;  " 

35 

* 

’Indicates  (a,  x)  not  in  "zone  ot  computation"  (see  Section  3). 
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3.  ZONK  OF  COMPUTATION 


The  function  P(a.  x)  (Q(a.  x>)  lias  the  property  that  for  a given  < > 0,  there 

exists  a function  x(a)  (x(a)>  such  that  I*  t (0  < f)  for  all  x ^ x(a)  (x  ^ x( a ) ). 

Of  course  x and  x are  also  functions  of  i f igure  _ .,nows  the  curves  x(a)  and  x(a)  on 
two  scales  for  t = c,  and  e = c, . Referring  to  the  figure,  I’  or  (J  is  computed  using 
( 10)  (28)  if  (a.  x)  is  contained  between  x and  x.  We  call  this  region  the  “zone  of 

computation."  If  (a.  x)  is  asserted  to  be  outside  this  region  by  the  estimates  given 
below  for  x and  x,  then  P is  set  to  zero  or  one  and  Q to  one  or  zero.  As  can  be  seen 
from  figure  2.  the  "zone  of  computation"  makes  up  only  a small  part  of  the  plane 
region  x i*  0,  a > 0.  Therefore,  it  is  advantageous  to  obtain  close  estimates  for  x and  x 
which  can  be  used  efficiently  by  the  program. 

Two  different  approaches  are  used  to  find  two  different  sets  ol  estimates  for  x 

and  x.  Associated  inequalities  are  given  which  can  be  easily  incorporated  into  the 

program  to  effectively  decide  if  (a.  x)  is  not  between  x and  x 


Let  the  first  estimates  for  x(a)  and  x(a)  be  denoted  by  x(  and  xr  respectively. 
The  expression  for  is  derived  in  [ 1 2 1 . It  is  given  by 


x,  = a 


(l  -1  + y v*5a)'.  i > 15 


where  y satisfies  the  equation 


erfety  \^2  ) = 2 ( 


The  constant  y is  determined  once  c is  specified,  l or  example,  when  c = e , y = 
7.1306  ■ • • . t - e,,  y = 4.892.  t = c^,  y = 3.29053.  A set  of  correction  factors  to 
improve  (32)  has  been  given  in  |8| , and  they  show  that  5^  > "x.  Therefore. 


Similarly 


Q(a,  x)  < r for  all  x > x > x . 


x ( = a 1 1 — y/\/9a  j , a > 15. 


where  again  can  be  established  from  |8)  that  x(  < x.  Therefore, 


P(a,  x)  < f for  all  x < x,  < x . 


In  order  to  use  these  results  in  the  program,  new  functions  I and  F are  defined 
Using  the  fact  that  Q is  a decreasing  function  of  x and  an  increasing  function  of  a,  we 
conclude  from  (32),  for  given  (a,  x).  that  if 


F = x 


a > 15 


(34) 


then  Q < t.  In  the  same  way,  we  conclude  from  (33).  that  if 


r = x 


• - r-  \ ' 

— y/v^a  ) ^ 0.  a * ) 5 
da  1 


(35) 


then  F < c.  Clearly  the  quantities  F and  F are  cheap  to  compute,  requiring  only  a 
square  root  and  a few  arithmetic  operations.  However,  they  are  only  useful  if  a is  not 
small  (a  > 15). 


The  other  estimates  for  x and  x arc  dei.  ed  by  x,  and  x,  respect  /c I y . We  first 
show  for  fixed  a and  c,  that  there  exists  a unique  positive  value  of  x.  say  x, . such  that 

R(a,  v,  ) = a ( I x,  (a  + I >)t . x,  < a.  ( Soc  <11*1.  (36) 

and 

l’(a.  x)  < t for  all  x x,  < a.  (37) 

Then  we  show  that  if  x,  is  determined  (uniquely)  from 
R(a,  xj  = x,  |(x;  + 2 a)/(x,  + I )|r.  0 < a s*  1.  or  a s 2 . (38) 


then 


Q(a,  x)  < c for  all  x J*  x,  . (39) 

For  1 < a < 2,  it  can  be  shown  Q(a.  x)  > e for  all  x x, . therefore,  we  replace 
(38) by 

R(a,  x2 ) = (x2  + I a ) c K a<  2.  x,  > a (40) 

and  show  that  (39)  holds. 
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Consequently,  since  £ is  an  increasing  function  of  x when  x < a,  . 37)  must  hold,  i.e., 

P(a,  x)  < $(x,  a)  < £U2 , a)  = e , x < x2  a.  (43) 

The  following  argument  for  the  existence  and  uniqueness  of  x2  in  (40)  depends  on 
the  function 

g(x)=R(a,x)  (x+1  a)e,  x>a,  a>l 

where  a and  e are  fixed.  The  proof  is  similar  to  the  one  for  (36).  Clearly  g(°°)  = —<*>, 
and  it  is  easy  to  show  g'(x)  < 0 for  x 3*  a.  The  proof  will  be  complete  if  it  can  be 
proved  that  g(a)  > 0.  We  have,  using  (26)  again  with  a i*  1,  that 

e » / • 2 a _ f 

The  expression  on  the  right  also  appeared  in  (41).  Therefore,  g(a)  > 0 provided 
(l/\/2ir  ) e '/|2  > e,  i.e.,  for  e < 0.367. 


g(a)  = R(a,  a)  e > 


We  proceed  to  show  that  (39)  holds  for  (40).  From  (2),  with  a - i < x, 


i r00 

Q(a.  x)  = — / exp  ( t + (a  1)  fn  t]  dt  , a > 1 (44) 

r(a)  A 


1 /-°°  1 + (a  - 1 )/t  , , 

< / exp  [ - t + (a  - 1 ) i n t]  dt 

F(a  )J\  - 1 + (a  l)/x 


e * xa  _ R(a,  x) 
f(a)(x  + 1 - a)  (x  + 1 - a) 


x > a - 1 > 0. 


It  is  easy  to  show,  by  differentiation,  that  R(a,  x)/(x  + 1 - a)  is  a decreasing  function 
of  x for  x > a - 1 > 0.  Hence, 


Q(a,  x)<  R(a,  x)/u+ I -aK  R(a,x2)/(x2  + 1 -a)  = e,  0<a<x2<x.  (45) 

Next,  we  establish  the  existence  and  uniqueness  of  x2  in  (38).  Consider  the 
function 


h(x)  = [(x  + 1 ) R/x(x  + 2 - a)j  - e,  0 < a < 1 , or  a > 2 and  x > (a  - 2), 


"J"  I""  ™ 


so  that 


h ( x ) = l R/x(x  + 2 


(x  + 1 ) (x  + 3 a)/(x  + 2 a) 


= (R/x) 


(I  a )( 2 a) 
+ 

x(x  + 2 a)2 


Now  if  0 < a < 1,  then  h(0  + ) = °°,  h(°°)  = e and  h(x)  < 0.  Hence,  there  exists  a 
unique  positive  value  of  x,  say  x,  such  that  h(x, ) = 0.  Clearly  this  implies  (38)  is 
satisfied  by  x,  since  x + 2 a > 0. 


If  a > 2,  x > a 2,  then  h(a  - 2 +)  = h(°° ) = c and  h (x)<  0.  Hence,  by  the 

same  arguments  as  above,  there  exists  a unique  positive  value  of  x,  say  x2  that  satisfies 
(38). 


For  0 < a < 1,  it  is  shown  in  (41  that  Q < h(x)  + e. 

This  is  easily  established  here  by  considering 

M(x)  2 Qfa,  x)  (h(x)  + c) 


and 


M'(x)  = ( 1 - a)(2  - a)  R(a,  x)/(x(x  + 2 a)] 2 

Then  if  0 < a < ! M(0  +)  = «>,  M(»)  = 0,  M'(x)  > 0.  Consequently  M(x)  < 0 for 

0 < a < 1 . Similarly  for  a > 2,  x > a - 2,  M(a  - 2 +)  = M(°°)  = 0,  M'(x)  > 0 so 

that  M(x)  < 0 for  a > 2,  x > a - 2 (note  also  that  M(x)  = 0 for  a = 1 and  a = 2). 

0 _ 

Finally,  since  h(x)  < o,  h(x)  + e is  a decreasing  function  of  x,  so  that  for  all  \>\2, 

Q(a,  x)  < (x  + 1)  R/[x(x  + 2 a)]  < (x2  + I)  R(x2 , a)/[x,(x?  + 2 a)]  = e (46) 

This  establishes  (39)  for  x2  satisfying  (38)  and  includes  the  values  a = 1,  a = 2.  Inci- 
dently,  h(x)  + e can  be  obtained  as  an  approximation  for  Q from  (98)  with  n = 1. 

In  (6,  p.  70],  the  first  inequality  in  (45)  is  derived  in  a different  way.  The 
inequalities  (34)  and  (35)  can  be  sharpened  [8] . The  same  can  be  said  for  (43),  (45), 
[3],  and  (46),  [4],  We  did  not  feel  that  sharper  estimates  at  the  expense  of  more 
computing  time  were  justified  since  the  above  estimates  were  already  very  good  as  can 
be  seen  by  looking  at  Table  2.  The  table  also  shows  that  (34)  and  (35)  are  better  than 
(43)  and  (46)  for  extremeiy  large  a. 


i 


i 


18 


10 
12 
18 
5.021 
100.018 
3 
2 


15.001 

25.001 
50.004 
75.006 

100.005 

250.007 

500.01 

1000.01 

5000.01 

10000.04 


1 

14.91 

15.019 

37.95 

24.92 

25.023 

33.53 

49.95 

50.038 

88.75 

74.97 

75.049 

121.57 

99.99 

100.055 

153.15 

250.05 

250.102 

331.68 

500.10 

500.147 

613.79 

1000.18 

1000.218 

1159.19 

5000.48 

5000.514 

5350.72 

10000.67 

10000.715 

10494.2 

14.73 

15.003 

50.49 

24.75 

25.004 

68.69 

49.78 

50.009 

108.65 

74.80 

75.015 

145.08 

99.82 

100.015 

179.69 

249.88 

250.029 

371.29 

499.95 

500.048 

668.09 

999.98 

1000.074 

1234.22 

5000.10 

5000.173 

5513.12 

10000.30 

10000.380 

10722.2 

4.93  (-7) 

4.94  (-7) 
4.92  (-7) 
4.91  (-7) 
4.89  (-7) 
4.88  (-7) 
4.85  (-7) 
4.85  (-7) 
4.85  (-7) 
4.83  (-7) 


4.21  (-7) 
4.43  (-7) 
4.66  (-7) 
4.75  (-7) 
4.82  (-7) 
4.91  (-7) 
4.95  (-7) 
4.97  (-7) 
4.99  (-7) 
4.99  (-7) 


5 (-7)  = e2 
5 (-7)  = e2 
5 (-7)  = e2 
5 (-7)  = e2 
5 (-7)  = e2 
5 (-7)  = e2 
5 (-7)  = e2 
5 (-7)  = e2 
5(-7)  = e2 
5 (-7)  = e. 


5 (-13)  = e 
5 (-13)  = e 
5 (-13)  = e 
5 ( 13)  = e 
5 (-13)  = e 
5 (-13)  = e 
5 (-13)  = e 
5 (-13)  = e 
5 (-13)  = e 
5 (-13) = e. 
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Table  2.  Comparison  of  x.(a),  x.(a)  with  x(a),  x(a),  i = 1,2.  (Continued) 

x,  x a Q(a,x2)  Q(a,  x, ) Q(a,  x) 


* 14.995 

24.992 

49.99 

74.975 

99.975 

249.975 

499.975 
999.97 

4999.99 

9999.97 


15 

25 

50 

74.994 

99.985 

249.99 

499.985 

999.98 

4999.99 

9999.99 


16.38 
25.58 
50.25 
75.15 
100.10 

249.99 
499  92 
999.83 
4999.54 
9999.32 


36.73 

51.67 

75.95 

100.68 

250.29 

500.15 

1000.06 

4999.92 

9999.83 


14.994 

24.989 

49.979 

74.957 

99.948 

249.906 

499.860 

999.781 

4999.497 

9999.277 


24.999 

49.999 


1.15 

4.99  (-7) 

1.27  (-7) 

5 (-7)  = e2 

5.36 

4.99  (-7) 

3.05  (-7) 

5 (—7)  = e2 

20.03 

4.96  (-7) 

4.21  (-7) 

5 (—7)  = c2 

37.17 

4.95  (-7) 

4.51  (-7) 

5 (-7)  = e2 

55.57 

4.94  (-7) 

4.66  (-7) 

5 (-7)  = e2 

177.02 

4.90  (-7) 

4.87  (-7) 

5 (-7)  = Cj 

394.90 

4.88  (-7) 

4.93  (-7) 

5 (-7)  = e2 

849.49 

4.85  (-7) 

4.96  (-7) 

5 (-7)  = e2 

4657.97 

4.83  ( -7) 

4.98  (-7) 

5 (-7)  = e2 

9514.47 

4.83  (-7) 

4.99  (-7) 

5 (-7)  = c2 

2.60  (-5) 

4.99  (-13) 

NA 

5 (-13)  = e, 

0.3075 

4.99  (-13) 

3.10  (-18) 

5 (-13)  = e, 

9.46 

4.99  (-13) 

1.23  (-13) 

5 (-13)  = e, 

22.85 

4.99  ( 13) 

1.39  (-13) 

5 (-13)  = e, 

38.15 

4.98  (-13) 

3.20  ( — 13) 

5 (-13)  = e, 

146.42 

4.96  (-13) 

4.36  ( 13) 

5 ( 13)  = e, 

349.58 

4.95  ( 13) 

4.70  ( — 13) 

5 ( 13>  = e, 

783.43 

4.94  ( 13) 

4.85  (-13) 

5(  13)  = e, 

4504.52 

4.92  ( 13) 

4.96  (-13) 

5 (-13)  = e, 

9295.57 

4.91  ( 13) 

4.97  ( 13) 

5 ( 13)  = e, 

*(34)  not  applicable  for  a = 2.60  ( 5)  = 2.60  x10  5. 
**x,  computed  from  (40)  instead  of  (38). 


As  stated  before,  the  use  of  (34),  (35),  (43),  (45)  and  (46)  contribute  significantly 
to  the  overall  efficiency  of  the  program.  Moreover,  the  computation  of  R(a,  x)  in 
(43),  (45)  or  (46)  is  not  always  wasted  if  these  inequalities  are  not  satisfied,  because 
this  function  is  also  needed  for  computing  P from  (10)  or  Q from  (12)  or  (15). 

In  the  next  section  we  give  a somewhat  detailed  discussion  of  Takenaga’s  method, 
for  use  when  a > 100,  and  the  changes  we  made  to  it. 
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4.  TAKENAGA’S  METHOD  AND  MODIFICATIONS 

In  this  section  Takenaga’s  method  [9]  for  evaluating  P when  a is  not  small  is 
summarized.  Changes  were  made  in  his  analysis  in  order  to  make  it  efficient  for  calcu- 
lation. In  addition,  it  was  extended  to  include  the  direct  computation  of  Q.  These 
modifications  will  be  given  in  detail.  Basically,  the  procedure  is  based  on  expansions 
in  terms  of  incomplete  normal  moment  functions,  [7] , and  it  works  best  when  a and  x 
are  close  and  large.  Even  though  the  method  can  be  used  for  a as  small  as  20,  it  is 
called  in  the  program  only  when  a > 100,  (and  3a  < 4x  < 5a).  Therefore,  it  will  be 
assumed  in  discussing  the  method  that  a > 100. 


t = |(z  + a)/0]  3.  (47) 

where  a and  (3  are  parameters  to  be  specified.  With  z as  a new  integration  variable  and 
a replaced  by  b + I,  (3)  transforms  to 


>(b  + 1 , x)  = 


Jl~ s 3 / Z + ai 
r 0 \ is  I 


(9)’ 


where 


r = - a . s = 0x,/J  - a,  b = a I. 


We  can  also  write  from  ( 1 ), 


P(a,  x)  = P(b  + 1 , x)  = 


-f, Sib- 


z)  dz  , 


where 


(3  /Z  + a\3b*; 
S(b.  z,  = l - — ) 


(z  + a\3  I 

(t)  f 


r(b  + l) 


± 


Then 


3 

log  S(b,  z)  ~ log  - 


P 


+ (3b  + 2)  log  - + (3b  + 2)  log  (1  + z/a) 

0 


(52) 


- |b  + - | log  b - b + - log  277  + L(b)  , 

where  the  quantity  in  brackets  of  the  last  line  represents  the  asymptotic  expansion  of 
log  l’(b  + 1 ) = log  b + log  T(b),  (see  (26)).  More  specifically, 

log  F(b  + 1 ) ~ |b  + - j log  b b + - log  2t7  + L(b)  , (53) 


with 

1 1 1 1 1 

L(b)  = + + ■ • • . (54) 

12b  360b3  1260b5  1680b7  1188b9 

(In  (27 7,  we  have  denoted  the  first  four  terms  of  L(a)  by  L(a)).  In  (52),  we  make  the 
substitution 


log(l  + z/a)  = - Y]  - ( - z/u)k,  1 < z/a  < 1.  (55) 

k = i k 

This  series  can  be  used  for  the  range  of  values  of  (a,  x)  for  which  the  method  is  to  be 
applied.  Indeed,  taking 


a = 3(b  + 2/3)' /2 

(56) 

p = 3(b  + 2/3),/6 

(57) 

and  noting  from  (47)  that  z is  an  increasing  function  of  t,  we  have  from  (49) 

1 < z/a  < (p /a)  x 1,3  - 1 = (x/(b  + 2/3)1 1/3  - 1 = s/a.  (58) 
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In  the  program,  we  use  (50)  only  for  cases  when  x < b + 2/3,  (a  corresponding  expres- 
sion for  Q(b  + 1 , x)  is  used  when  x > b + 2/3  and  will  be  discussed  later  in  this  section). 
Consequently,  for  x < b + 2/3,  (58)  shows  z/a  < 0 so  that  the  right  hand  inequality  in 
(55)  is  satisfied,  (it  is  interesting  to  note  that,  regardless  of  the  constraint  on  x, 
z/a  < 1 provided  0 < x < 8(b  + 2/3)  which  always  holds  for  a > 6 and  e < P < 1 - e, 
where  e,  < e . This  follows  from  the  results  of  the  previous  section). 

The  left  inequality  in  (55),  on  the  other  hand,  cannot  be  satisfied,  because  z takes 
the  value  - a when  t = 0.  This  difficulty  is  circumvented  in  the  following  way.  Let 


i,  x)  = P(a,  X)  + J e-'  tb  dt/r (a)  , 


where  a = b + 1,  as  noted  above,  and 


X-a  /l  - -y/3v^rj3,  y = 1 1 .0  . (60) 

Then  by  the  results  of  the  previous  section  (see  (33)  with  y replaced  by  y)  P(a,  X)  < 
4 x 10  28 . Hence,  we  can  take  the  second  term  on  the  right  hand  side  of  (59)  in  place 
of  P(a,  x)  with  an  insignificant  loss  in  relative  accuracy.  Under  these  conditions,  we  get 
from  (47),  that  when  t = X,  z = r,  and 

r !a1/3  |l  — - y/3\/a  j - a (61 ) 

= a J (a/(a  1/3)1 1/3  ^ - v/3v^)  1 j 

= a|(l+  — + — -L+  -j(l  1 y/3  Va' ) I | 

| \ 9a  81  a2  M 9a  ' I 

>aW|+l  + _L]/l  _L  y/3  \/a  ) l[>  -^(1+-  + ^;) 

I ' 9a  81a2'  ' 9a  ' I 3^  ' a 81a2' 

Thus,  since  in  our  case  a > 100,  we  see  from  (61)  and  (58)  that 

— - — (l  + — + — — ) < r/a  < z/a  < s/a  < 0. 

3 v^a  ' 9a  81a2'  " 


This  justifies  the  use  of  (55),  since  for  y = 1 1.0  and  a > 100 


y 

3\/a 


1 

1 + — + 
9a 


< 0.37 


The  choices  for  a and  (3,  given  by  (56)  and  (57),  were  shown  by  Takenaga, 
after  using  (55)  in  (52),  to  result  in  zero  coefficients  for  the  terms  z/a,  (z/a)3,  logb, 
b in  (52).  Thus,  with 

log  a/0  = (1/3)  log(b  + 2/3)  = (1/3)  [log  b + log(  1 + 2/3b)l , 


he  found 


z2  a2 


log  S(b,  z)  = log  C - - log  2?r — y\  ( 1/k)  (-  z/a)k , 


k ■ 4 


(62) 


where 


log  C = (b  + 1/2)  log(l  + 2/3b)  2/3  L(b) 


(63) 


= (b+  1/2)  £ (l/k)<  2/3b)k  2/3  - L(b) 


111  2 
+ + 


691  16 

+ 


36b  81b2  360b J 1215b4  306180b5  15309b6 


373  80 

+ 


44069 


3674160b'  177147b8  38972340b9 


Then  from  these  results,  he  obtained 


S(b.  z)  = 0 


V5T 


e_/  12  exp 


— V ( 1/k)  ( z/a)k 

3 ^ 


k 4 


(64) 


*The  coefficients  in  (63)  and  (65).  or  (68)  and  (80),  were  computed  by  Alfred  Morris 
with  his  “Flap”  Algebraic  Routine,  [5] . Additional  coefficients  are  given  in  Appendix 
B.  In  ( 9 1 there  is  a misprint  in  (14)  and  another  in  (19)  of  that  paper. 
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Now  expanding  the  second  exponential  function  of  (64)  in  powers  of  z/a,  one  gets  the 
power  series  in  (65),  namely 


2 00  J 1 / 1 \ 

exP  - — 52  (l/k)(  z/a)k  =1  z4  (z6 z8) 

3 12a2  18a4  \ 16  ' 


1 

/z»  1L  z'o 

+ — z12  \ 

j 

24a6 

\ 225 

432  / 

1 

1 

I24 

'78 

2.3414  67226  5 1 062  x 

10  22  z48) 

4. 

1 

/5+_L  lZ- 

+ 1 L* 

29 

7i  1 

+ iz.) 

15a3 

3a5  ' 7 

60  / 27a7  \ 

140 

160  / 

4. 

1 / 

1 zn  193 

.,  1187 

-,13+  2' 5 

1 

z.?\ 

T 

j 

a9  V 

33  22680 

2268000 

155520 

Z ) 

) 

» 

+ 

1 

(-z27 

'81 

2.2478  08537  45020 

x 10  21 

z49)  + • • 

In  order  now  to  carry  out  the  integration  of  (64),  using  (65),  we  need  to  evaluate 
integrals  of  the  form 


H 


B = I z2'  Z(z)  dz,  A = 


■JT-~ 


Z(z)  dz. 


where  j is  a non-negative  integer  (in  the  program  j < 24)  and 

1 2 

Z(z)  = — e'1  12  . 


In  [9] , the  indefinite  integration  of  (64),  using  (65),  is  carried  out,  and  a cumbersome 
procedure  is  suggested  for  evaluating  the  constant  of  integration.  The  approach  below 
is  more  direct  and  easier  to  apply. 
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Wc  express  the  integration  of  (64)  in  inverse  powers  of  a. 


s 

Sib.  /.)  dz  = C 


« 


k 


(6") 


where 


d 


o 


I _ 


1 - 


and  in  general  for  k > 2. 


Ik  even.  N()u  I) 


Ik  odd) . 


(68) 


The  eonstants  N'kj,  Mkj  which  result  from  the  expansion  in  *65)  are  given  in  Appendix 
B.  Some  particular  values  are  N,,  = -p.  M?  5 = 24>/( 27  x 140).  If  k is  even, 

a typical  term  of  the  sum  on  the  right  hand  side  of  (68).  multiplied  by  « " k.  is  given  by 


a k N.  B h N.  u k f r>  Z dz  N.  u k fL  z:j 

k|  j k)  I k i I 

-/  JT  J OO 


Z dz 


(69) 


or  in  the  case  of  k odd. 


u k M.  A - M 

kj  I k| 


u k f r'"  Z dz  Mkj  a k fL  z2^'  Z 

J <.TO  J OO 


dz 


(70) 


We  want  to  show  that  the  last  term  on  the  right  hand  side  of  these  identities  can  be 
dropped  for  b(=  a I ) > 99  with  no  significant  loss  in  relative  accuracy.  First,  we 
obtain  an  upper  bound  on  r.  From  (61). 


r = a 


| / b + 1 \ '/3 
Mb  + 2/3  ' 


9,b  + 1 ) 


b + 1 
b + 2/3 


3x/F+i 


I 


Then 


1 b + i y/J 

1 

| 

1 

1 + 

1/3 

1 

1 

'b  + 2/3  ^ 

9(b  + 1) 

3(b  + 2/3) 

9(b  + 2/3)[  1 + l/3(b  + 2/3)] 

■H) 


3 \‘/3 


a2  ( 1 + 3/a2) 


(71) 


< 


/ 1 1 5 W 13  9 27 \ 

1 +-  - + — 1 — - + — 

\ a2  a4  3a6/  \ a2  a4  a6  a8/ 


3 9 27' 


1 10  1 40  1 

<1+  — — — + — — <1  + 1 /a4 , (a  > 1 00) . 

a4  3 a6  3 a8 


Therefore. 


_ ib  + J/Jv'/o  1 

r < y +-<  10.9938  (b  3*  99,  y=11.0). 

V b + 1 ' a3 

Now,  letting  v = 10.9938  and  taking  k = j = 0 in  (69),  we  have 

f l Zdz  = f Zdz<  rZ  dz  < 2.05  X 10  28  (See  [l.p.  972)) 

J —oo  J r -'v 

Consequently,  the  last  term  in  (69),  with  k = j = 0,  can  be  neglected,  (the  constant  C, 
in  (67),  does  not  affect  this  conclusion,  since  it  is  always  close  to  one,  for  a > 100.  as 
(18)  or  (63)  shows). 

New  consider  the  absolute  value  of  a particular  term  in  (67)  for  k > 0 and  even. 
|NkJ  a k z:>  Zdz<  |Nk  | a'k  f °°z2j  Z dz . 

oo  J v 


Then  with  u = z2/2  and  the  use  of  (45), 


28 


TWC "XIK1 1 1 J 


I unp^  I,  iiiu,ii>»jih,>».  iptiy;! 


|Nt,|  « 1 f“r'  Z .12  = |\,|  ^Ll  f“  u-  e-  du 


(72) 


= |NJ  T~?=-  QG  + 1/2.  v2/2)  rG  + 1/2) 


- \ n 


|\  | u k 21  e v‘/2  (v2/2)J+ 1/2 


2\/j? 


l(v2/2 ) + 1 j 1/2] 

l\.|«  k v^* 


(See (45)) 


v-  VI 


|v:  + 1 2 j ] 


A similar  analysis  lor  k odd  shows 


r“  a k e v /:  v2jt2 

|M  | u k / z2)M  Zd/.  *|MJ  — ■ 

1 k)l  Jv  1 kjl  y/Tn  (V2  2j) 


(73) 


It  can  be  shown  lor  every  case  considered  in  the  program,  i.e..  0 < k < 24, 
(k  + 2)/2  < j < k (k  even),  (k  + I )/2  < j < (k  1 )(k  odd)  that  by  bounding  the  right 
hand  sides  of  (72)  and  (73).  the  left  hand  sides  of  these  inequalities  are  very  small. 
For  example  setting  b = 99  (for  winch  u k takes  its  largest  value),  k = 23,  j = 22. 
M,3  ,2  ^ 3 \ 10  19  (see  Appendix  B).  we  have  from  (73) 


M 


a 23  f 
:j.::  / 


> Z d,  < Ill'Ll-  » 10  M x 4.5  10”  - OO0-”). 
\/2rT  x 76.86 


For  the  same  value  of  k with  j = 12.  M, , ' 0.0133 


M 


23,12 


Jr-  oo 

z25 

V 


I A \ 1 0 2 x I . I x 1 0 34 

/ dz  < x 6.7  = 0(10  37). 

x 95.86 


Thus,  the  second  tenns  on  the  right  hand  sides  of  (69)  and  (70)  will  be  discarded, 
or  equivalently  the  integrals  Bj  and  Aj  in  (66)  are  replaced  by 


* r 


B,  = / z2>  Z dz 


s-/>- 


Z uz 


(74) 


(75) 
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These  integrals  are  easily  evaluated  on  a computer  by  using  the  recurrence  rela- 


Bj  = (2j  1 ) Bj_ , s2^- * Z(s) 


A . = 2j  , s2j  Z(s) 


/~s  1 

= / Z dz  = - crfc  (Isl/VT  ),  s < 0 
J - OC  2 


.-*12  = 


There  is  no  problem  to  compute  BQ  since  a very  efficient  and  accurate  error  func- 
tion routine  is  available.  It  is  based  ori  Cody’s  minimax  approximations,  (2),  and  it 
yields  14  significant  digits  for  erf  or  erfc  depending  on  the  magnitude  of  the  argument. 
For  completeness,  the  particular  Cody  approximations  that  we  use  are  specified  in 
Appendix  A.  A discussion  of  how  s is  computed  is  given  on  page  41. 


Finally  the  quantity  C,  whose  logarithm  is  given  by  (63),  is  expressed  as  a power 
series  in  inverse  powers  of  b in  order  to  evade  computing  log  C from  the  series  in  (63). 
and  subsequently  C from  exp  (log  C).  Thus  we  obtain 


C = 1 + — b 
36 


31  , 3413 

b 2 + b 

2592  1399680 


361733 

+ b 4 

201553920 


113888281  c 7565202533 

b 5 + b 6 

50791587840  7836416409600 

Takenaga's  procedure  is  extended  to  advantage  by  using  the  same  approach  when 
x > b + 2/3  = a 1 /3  to  compute  Q instead  of  P.  As  noted  earlier.  1’  in  this  case  can  be 
close  !o  one  so  that  0 must  be  computed  directly  to  maintain  the  specified  relative 
accuracy. 


*See  footnote  on  page  25. 
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The  computation  of  Q only  requires  changing  the  limits  of  integration  in  (74)  and 
(75)  from  - «>  and  s to  s and  «,  respectively.  Thus  (76)  - (79)  are  replaced  by 


B.  = (2j  1)  B.  , + s2j-‘  Z(s) 


\ = -)  A.  , + s2j  Z(s) 


B = — — f e ,2/2  dt  - - erfc  (s/\/T) 

° V2tJs  2 


A0  = ~ e's2/2  = Z(s) 
V 2n 


(81) 

(82) 

(83) 

(84) 


Certainly  these  changes  are  quite  simple.  We  also  note  from  (49)  or  (58)  that  s > 0 
since  x > b + 2/3. 

In  order  to  establish  these  results  for  Q(b  + 1,  x),  we  require  the  integration  of 
S(b,  z)  in  (51)  from  s to  with  s > 0.  Again  this  raises  the  question,  when  we  pro 
ceed  as  before,  concerning  the  validity  of  using  the  expansion  in  (55).  In  order  to 
justify  its  use  we  must  show  that  |z/a|  < 1 for  all  z in  the  domain  of  integration. 
Clearly  this  will  not  be  the  case  if  z e [s,  °°] . We  get  around  this  problem  much  like 
before.  Note  that  since  x > b + 2/3,  a lower  bound  on  z/a  is  easily  found,  namely 

(z/a)  > s/a  = (x/(b  + 2/3)]  1/3  1>0 

An  upper  bound  on  z/u  is  found  by  considering 

X 


Q(a,  x)  = Q(a,  X)  + 


X" 


r dt  , a = b + 1 , 


(85) 


where  X = a(  1 - — + y/3\/a  )}  > x > 0.  It  follows  from  (46)  and  [ 1,  p.  972)  that  for 

v = 1 1 .0,  Q(a,  X)  < 7 x 10  2‘'.  Thus,  we  can  replace  Q(a,  x)  by  the  second  integral  on 
ihe  right  hand  side  of  (85).  We  get  an  upper  bound  on  z/a  by  using  t < X in  (47). 
Indeed, 


z < 


a |(/j/a)X*/->  1]  =a  (l  - - +7/3^)  - 1 

Mb + 2/3/  \ 9a  ! \ 


/b  + 2/3  \ 1 ^ 

<( y + l/aJ  < 10.9939  = v,  (See  (71)). 

V b + 1 f 
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Since  b > 99,  we  have,  by  noticing  that  v/a  is  a decreasing  function  of  b, 


0 < s/a  < z/a  < v/a  < 0.37. 


This  result  allows  us  to  use  (55)  when  s > 0,  (x  > b + 2/3). 


The  algebraic  procedures  vhich  yield  (65)  in  inverse  powers  of  a,  and  the  subse- 
quent integration  show  that  ino . -'idual  terms  contain  integrals  of  the  form  (after 
dropping  Q(a,  X)  in  (85)) 


B = 


r 


Z dz, 


Z dz. 


(86) 


where 


Z 


\/2n 


-S 12 


Proceeding  as  in  the  case  for  x < b + 2/3,  the  expressions  corresponding  to  (69)  and 
(70)  are  given  by 


a k N.  B.  = N 

kj  j kj 


Xoo  jo 

z2j  Z dz  - Nkj  a " k J z2i  Z 


dz 


(87) 


Jp  OC  OC 

z2>"  Z dz  M a k / z2j+1  Z 

s kj  Jv 


dz . 


(88) 


We  note  thu*  the  coefficients  of  (68),  Nkj  and  Mkj,  are  unchanged  and  that  (68) 
holds  wit!  B^  and  A;  given  by  (86).  But,  we  have  argued  above  (see  (72),  (73))  that 
the  second  term  on  the  right  of  (87)  and  (88)  can  be  dropped  or,  equivalently,  the 
integrals  Bf  and  /V  in  (86)  can  be  replaced  by 


z2’  Z dz 


Jr-  oo 


Z dz. 


(89) 


It  is  interesting  to  note  tha*  these  integrals  were  obtained  directly  by  formally  integrat- 
ing (64)  from  s to  00 . They  are  evaluated  by  the  recurrence  relations  (81)  and  (82). 


The  final  integrated  expressions  for  P(b  + I , x)  and  Q(b  + 1 , x)  are  given  by  ( 1 7) 
and  again  in  (115)  ot  the  next  section  where  the  computer  program  is  discussed. 
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5.  DFSCRIPTION  OF  PROCR  \M 

In  this  section  we  summarize  the  equations  and  inequalities  on  which  the  C DC- 
6700  computing  program  is  based.  In  addition,  the  general  features  of  the  program  are 
described  and  a master  flowchart  is  included  at  the  end  of  this  section  as  Figure  3. 

Ihe  program  requires  as  input  5 quantities,  the  two  locations  where  two  argu- 
ments a and  x are  stored,  the  two  output  locations  for  P and  0.  and  a location  which 
is  used  to  specify  whether  c,  (=  5 x 10  1 3 ),  c2(=  5 x 10' 7 ),  or  e3(=  5 x 10  4 ) is  to 
be  used. 


If  c ( c _,  ) If.-  1 is  specified,  then  the  output  quantities  P and  Q are  given  correctly 
to  within  one  unit  in  the  twelfth  (sixth)  i third!  significant  digit.  Of  course,  if  (a,  x)  is 
outside  rhe  "zone  of  computation  ' I*  and  Q may  have  no  correct  significant  digits 
since  tin  y are  simply  given  as  I and  0 or  0 and  1 . The  absolute  error  in  such  cases,  of 
course,  will  always  be  within  the  specified  c. 


In  order  to  give  some  estimate  of  the  average  computing  time  per  case,  two  grids 
of  points  in  the  ax-plane  were  used  for  c = t f and  two  others  for  c = e . The  two  grids 
for  e,  are  specified  first,  identifying  them  as  A,  and  Br 


(■rid  A, 


Cirid  B. 


a - 0.0001  (0.001)  0.05 
x = 0 (0.1 ) 126.1 


a =0.0001  (0.001)0.05 
x = 0 (0.1 ) X (a) 


a =0.1  (0.1)  3! 
x = 0(0.5)  1260 


a = 0.1  (0.1  ) 31 
x = 0(0.5)  X (a) 


a = 31.0  (0.5)09.5 
x = 0(0.5)  1260 


a = 31 .0  (0.5)  09.5 
x = x,  (0.5)  X (a) 


a = 100.1  (5)  1000. 1 
x = 0(10)  1260 

X (a)  minimum  grid  value  of  x.  given  a.  such  that  Q 


a = 100.1  (5)  1000.1 
x = Xj  ( 10)  X (a) 


x(  value  of  x.  given  a.  for  which  P(a.  x)  ~ t.  See  (33). 
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= a 


(l  1 y /3Va")3  , y = 7.1306,  (y  = 4.892,  e = e2 ). 


All  the  grid  points,  (a,  x),  specified  are  considered. 


Denoting  the  corresponding  grids  for  e = e 2 by  A,  and  B2 , we  have 


Grid  A2  e 

a = 0.0001  (0.001)  0.05 
x = 0 (0.1)  1 17.1 

a =0.1  (0.1)  20 
x = 0 (0.5)  1170 

a = 20.5  (0.5)99.5 
x = 0(0.5)  1170 

a = 100.1  (5)  1000.1 
x = 0 ( 10) 1170 


e.  Grid  B2 


a =0.0001  (0.001)0.05 
x = 0(0.1  )X(a) 

a = 0.1  (0.1)20 
x = 0 (0.5)  X (a) 

a = 20.5  (0.5)99.5 
x = x,  (0.5)  X (a) 

a = 100.1  (5)  1000.1 
x = x,  (10)  X (a) 


The  major  portion  of  the  points,  (a,  x),  in  the  A grids  lie  outride  the  “zone  of 
computation,"  whereas  almost  all  of  the  points  in  the  B grids  lie  within  the  “zone  of 
computation.”  The  A grids  were  designed  to  get  a representative  average  time  per  case 
in  the  situation  where  any  point  in  the  region  0.0001  < a < 1000.1  0 < x < U*  (where 
U denotes  the  maximum  value  on  x)  is  equally  likely  to  occur  as  input.  On  the  other 
hand,  the  B grids  were  set  up  to  obtain  a represen:  itive  average  time  per  case  in  those 
situations  where  the  input  | oint  (a,  x)  is  most  likely  to  occur  in  the  “zone  of  compu- 
tation.” In  Tabic  3,  the  results  for  the  various  grids  are  summarized. 


•We  note  for  a = 0.0001  (0.001)  0.05  U is  taken  as  126.1  for  e = e,  and  117.1  for 
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Table  3.  Summary  of  Grid  Results 


\FQ 

grid\ 

GO) 

(12) 

(13) 

(15) 

(16) 

(17) 

li 

TOTAL 

CASKS 

AVG. 

TIMH 

(msec) 

GRID  A, 

20021 

39089 

1891 

17520 

769 

5872 

1 131595 

1216757 

0.75 

GRID  B, 

19958 

39155 

1879 

17512 

769 

5860 

2274 

87407 

2.02 

GRID  A2 

12920 

20708 

561 

5606 

769 

4020 

876965 

921549 

0.68 

GRID  B2 

12788 

20783 

551 

5598 

769 

3993 

1687 

46169 

1.59 

The  numbers  in  parentheses  refer  to  the  equations  in  the  text,  e.g.,  (10)  was  used 
20021  times  in  Grid  A( . The  column  headed  K gives  the  number  of  cases  that  were 
recognized  by  the  program  to  be  outside  the  “zone  of  computation.”  The  average 
time  per  case  is  given  in  the  last  column  in  milliseconds.  The  much  smaller  computing 
times  for  the  A grids  are  to  be  expected,  since  most  of  the  cases  were  outside  the  “zone 
of  computation.” 


The  flow  of  the  program,  depending  on  the  input,  should  be  easy  to  follow  from 
the  flow  chart  at  the  end  of  this  section.  The  circled  numbers  are  used  to  denote  sub- 
programs for  computing  P or  else  Q,  and  the  number  itself  corresponds  to  the  equation 
in  Section  2 upon  which  the  subprogram  is  based.  For  example,  (K))  refers  to  that 
subprogram  which  uses  (10)  for  the  computation  of  P;  (l3)  refers  to  that  subprogram 
which  uses  (12)  to  determine  Q.  The  quantities  F and  F are  defined  by  (34)  and  (35). 
The  term  R(a,  x)  = e~  * xa/l'(a),  (see  also  (II)),  and  e refers  to  either  e, , e2 , or  e , 
(see  (6)  and  (7)). 


We  continue  with  a brief  description  of  the  algorithm  used  with  each  subprogram. 
It  is  assumed  (a,  x)  is  in  the  “zone  of  computation”  and  that  all  the  subprograms 
except  (n)  require  a < 100.  * 


If  a > \,  then  (l0)  is  used  if  a > Also,  (To)  is  used  when  2a  is  not  an  integer 
provided  1 < a<  x < fn  10ora<  l,x<  1.5  and  P < 0.90,  or  a < 1 , x ^ 1.5  and 


"If  a > 100,  (To)  is  used  if  3a  > 4x  and  (Q)  is  used  if  4x  > 5a. 


\ 
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R > 0.101  x (x  + 2 - a )/(x  + 1).  Let  sn|  denote  the  (n  - 1)  st  term  of  the  series 
in  (10),  then  the  n,h  term  is  found  from 

sn=Hr)sn-r  so  = l>  n = 1 > 2,  • • ■ (90) 

\a  + n/ 

The  summation  of  the  series  is  terminated  at  n = N when 
Sjj  < e/2  and  x/(a  + N + 1)  < 2/3, 


which  assures  a truncation  error,  ^ sk  < e,  (see  Appendix  D). 

N + I 

Program  (£j)  applies  when  1 < a < x,  and  \>  31  if  e = e, ; x > 17ife=e2; 
x S*  9.7  if  e = e3 . The  (n  + 1)  st  term  of  the  asymptotic  series  (12)  for  Q is  given  by 


! = (!_ V , V = 1,  n = 0,  1 , 

n + 1 y I n ’ 0 ’ ’’ 


The  summation  is  terminated  when 


I v„  + 1 1 < £ . (92) 

Program  (0)  is  used  when  a < x < 31  for  e = e, , a ^ x < 1 7 for  e = e2,  or 
a < x < 9.7  for  e = e3 , and  if  2a  is  an  integer.  The  algorithm  for  computing  Q from 
(13)  then  is  given  by 


W(k  + 


g l,x)=  ( )w(k  + g 2,  x) 

\k+g  1' 


Q(k  + g,  x)  = Q(k  + g l,x)  + W(k+g  - 1 , x),  k = 1 , 2,  ■ , a - g,  (94) 
where  g = 1 if  a is  an  integer  or  g = 1/2  if  a is  not  an  integer.  Also 

W(k  + g - 1 , x)  = R(k  + g 1,  x)/(k  + g - I ) (95) 


x)  1/2  e * g = 1/2 


W(g  - 1 , x)  s 


erfc  (Vx)  g = 1/1 


Q(g,  x)  = 


(See  (14)). 


> mi  wuij..i|jai»  "M,  IP*  no u " u 1 


Program  (0)  is  based  on  a continued  fraction  expansion  for  Q.  It  is  used  when 
2a  is  not  an  integer  and  1 < a < x with  £n  10  < x < 31  for  e = e ( (£n  10  < x < 1 7 for 
e = , or  fn  10  < x < 9.7  for  e = e3 ) or  a < 1,  x > 1.5  and  R «;  0.101  x(2  + x - a)/ 

(x  + 1).  The  algorithm  for  computing  two  successive  approximates,  based  on  (15),  in 
each  iteration  follows.  Let 


D,  = D = 1,  E = x,  E = (x  + 1 - a),  n = 1,  2, 


Then 


D.  , = x D,  +nD,  , 

2n»l  2 n 2 n - I 


(98) 


E.  , - x E.  +nE. 

2 n * I 2n  2n - I 


Increase  n to  n + 1 


D,_  = D2n_,  + (n  - a)  D 


2 n 


2 n - 2 


(99) 


F =F  + (n  ;> ) F 

L2n  L 2 n - 1 Ul  ;L2n-2 


The  procedure  is  stopped  when 


D^„  »2n  , 


E,  E, 

2n  2 n I 


< e 

P2n 

E2n 

(100) 


Program  is  used  when  a < 1 . x < 1.5  and  P > 0.90.  From  ( 1 6),  three  series 

are  evaluated,  J,  L,  H. 


We  have 


J = -a  Z : 


( x)k 


k x | 


(a  + k)  k! 


- - a£  Jk  = - a Z V(a  + k> 


(101) 


then 


Tk  = (-  x/k) Tk_ , , Jk=Tk/(a  + k)  k = 2,  • • • , 


T,  = - x,  J,  = 


a + 1 


(102) 

(103) 
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The  iterations  are  stopped  when 


'V<-«£  j, 


i I 


For  L,  we  have 


L = a fn 


1 + £ (a  fn  x)k/(k  + 1)! 


k - I 


= a i n x 


1 + E 


Then 


L> “ (ttt)  L‘  '•  L|  *<a,l”‘W2-  k = :-' 


These  iterations  are  terminated  when 

-Lki  < e/2 


The  H series. 


17 


IX  a 


0 < a < 1/2 


H = 


(1/a) 


1 7 


(1  a)  + £ Ck  (a  l)k-‘  1/2  < a < 1, 


(104) 


(105) 


(106) 


(107) 


(108) 


* 


is  terminated  when  k = 17  or  before  if  Kn  < (e/2) 


EK, 


, where  K.  = C ak  1 

’ k k 


if  0 < a < 1/2  or  Kk  = Ck(a  1 )k  ~ 1 if  1 12  < a < 1 . We  note  also  that  the  11- 

series  make  up  part  of  the  series  that  are  used  for  evaluating  l/r(X)  when  0 < X < 1, 
(see  (25)).  Consequently,  when  0 < a < 1,  11  is  obtained  during  the  computation  of 
1/C(a)  and  stored.  The  evaluation  of  1/P(a)  is  taken  up  later  in  this  section  on  page  42 
and  in  Appendix  A. 
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Program  © is  the  only  one  used  when  a > 100,  and  3a  < 4x  ^ 5a.  The  indi- 
vidual integrals  Ak , Bk  that  appear  in  (17)  are  given  by  (20)  and  (21)  if  x < a 1/3 
and  by  (23)  and  (24)  if  x > a 1/3.  In  case  the  first  inequality  is  satisfied,  P is  com- 
puted, and  the  recurrence  relations  for  evaluating  the  integrals  of  (21)  and  (20)  are 
given  by  (76)  and  (77).  In  case  the  second  inequality  holds,  Q is  computed,  and  the 
recurrence  relations  for  evaluating  (24)  and  (23)  are  given  by  (81)  and  (82). 

We  write  again  the  basic  equations  for  © , with  (76),  (77),  (81 ),  (82)  slightly 
changed  for  greater  efficiency  in  computation.  Let 


B = Z(s)  b A = Z(s)  a (z(s)  = — c”'*12 ) • 
1 ’ J 1 \ / 


(109) 


Then 


v2|- 


| b,  = (2j  Db  , s2 
1^  = 2]  a , s2) 

replace  (76)  and  (77)  and  are  used  wiien  x <a  1/3  with 

J bQ  = (l/2Z(s)|  erfcdsl/v?) 

If,  on  the  other  hand,  x > a 1/3,  then  using  (109),  (81)  and  (82)  become 

j b.  = (2j  1 ) b , *•  sJJ- 1 

| af  = 2j  , + s2j 


(110) 


(111) 


(112) 


with 


b0  = | l/2Z(s)]  erfc  (s/s/5) 


(113) 


ao  = 1 ' 


In  all  cases 


s = 3(a  l/3)l/2  j [x/(a  l/3)j'/3  l}.  (Sec  (22)). 


(114) 
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Thus,  in  terms  of  ak,  bk,  instead  of  Ak,  Bk,  ( 17)  becomes 


T(a,  x)  = Z(s)  C 


I I / 1 v 

b b,  b — b.) 

12a2  18a4  \ 3 16  7 


(115) 


I 


24a6 

1 


37  I 

b„  — b + — b 


225  s 432  h 


743  49 

h h + b. 


30a8  ' 5 3024  6 4320 


-6.) 

944  8/ 


82944 


1 / 1 

— — b,, 

a24  '78 


2.3414  67226  51062  x 10 


” b><) 


1 5aJ 


/a 

a4  \ 1 / 

29  1 v 

M 

ac  + a,_ 

\ 7 

60 1 27a  ' 

140  5 160  / 

1 

i a< 

193  1187 

1 V 

+ 

— 

a/  + a. 

ax  1 

a 

'33 

22680  2268000 

155520  ' 

--  (iu 

a28  '81 


M78  08537  45020  x 10  21  a,J 

2 4 / 


where,  with  b = a 


1 , 31  3413  , 361733 

('  = I + — b 1 b * + b ’ + b 4 (116) 

36  2592  1399680  201553920 


113888281  £ 7565202533 

b 8 + b " 


50791587840 


(See  (18)). 


7836416409600 


The  program  uses  no  more  terms  than  are  shown  for  (115)  and  ( 1 16).  The  numerical 
coefficients  in  these  equations,  and  some  additional  nes  for  (I  16).  not  used  in  die 
program,  are  given  in  Appendix  B. 
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The  right  hand  side  ol  ( 1 15)  is  computed  by  the  following  algorithm 

(a)  Compute  aQ,  b0,  [Q  (=  b0 ) (use  (111)  for  s «£  0 or  (1 13)  for  s > 0) 

(b)  Compute  1/a2,  1 la  (a  = 3(a  l/3)l/2) 

(c)  Set  i = 2,  - = <o  <-  denotes  accumulated  sum) 

(d)  Compute  b(  ),bj,ai  , , a(  (by  (110)  for  s < 0 or  (112)  for  s > 0) 

(e)  Compute  f = ( S \b,)U-  V-  = (E  M,+ . ,i  0 A'* ' 

'1+1/2  ''  ' I +i/2  " 

(See  (67)  and  (68).t 

(f)  Add  and  to  accumulated  sum,  I,  to  obtain 


II 


(g)  Is  l/j  < t b0,  and  1 ^ , i | ^ < b0,  or  is  i = 24? 

(It)  If  no  in  both  cases  then  increase  i by  2 and  return  to  (d) 

(i)  If  yes  in  either  case,  for  (g).  proceed  to  (j ) 

(j)  Compute  Z(s)  Cl  and  exit. 

The  calculation  of  s does  not  result  in  any  significant  increase  in  relative  error  in 
computing  I*  or  Q.  A simple  analysis  shows  the  relative  error  in  P or  Q is  not  greater 
than  twice  the  absolute  error  in  v.  Nevertheless,  in  order  to  minimize  the  absolute 
error  in  s.  it  is  evaluated  by  the  algebraic  equivalent  of(l  14),  namely 

».  /.<i  + 1— r + f-M"*!  ■ 

(a  13)/  \ \a  1/3/  'a  13/  ) 


where 


b = (x  a)  + 1 /3  . 


r 


1 


If  (10),  (12),  (15),  or  (16)  is  used  the  function  R(a,  \),  as  defined  in  (11)  is 
computed.  If  a > 30  its  computation  is  carried  out  by  using  (28)  with  L(a)  given  by 
(27).  This  eliminates  any  scaling  problem  with  evaluating  R(a,  x)  for  30  < a. 

If  a < 30,  then  R is  evaluated  by 

R(a,  x)  = [e  x+a  r"*|/r(a).  (117) 

where  the  following  procedure  is  used  to  compute  l'(a).  Let  X be  the  fractional  part 
of  a. 


If  X = 0,  then 

( I’(j  + i ) = j l'(j),  j = 1 , 2,  • • • , a - 1 , a>l. 

< (118) 
( I'd)  = 1.  a = I . 

If  X # 0,  then  1/I'(X)  is  found  to  14  significant  digits  by  the  polynomials 


i/m)  = 


I 7 


1 + 


Ec. 


* E c»<»-  '*1"' 


2 


0 < X < 1/2 


1/2  < X < 1 ; 


(119) 


the  Ck’s  are  given  in  Appendix  A with  the  rule  for  terminating  these  series.  Then 
f(a)  is  computed  from 


jr(j+X)  = (j  I + X)  r(j  - 1 + X),  j = 1 , • • • . a-X,  a>l 
I r(a)  = r(X),  a < 1 . 

For  values  of  a > 20  it  was  observed  that  the  asymptotic  series  for  l n r(a),  (26), 
yielded  a faster  algorithm,  in  spite  of  evaluating  P(a)  from  e^n  than  the  one 

based  on  ( 1 1 8)  - (1 20).  However,  the  relative  error  was  larger.  In  the  procedur : above 
1.5  more  significant  figures  were  obtained  in  some  cases,  for  c = e, , than  by  using  the 
asymptotic  series. 
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The  error  function,  erf(x),  and  its  complement  e»fc(x)  are  evaluated  by  Cody’s 
minimax  approximations,  [2] , (see  Appendix  A).  We  note  that  the  error  function,  or 
its  complement,  is  needed  in  ( 1 T),  (111),  (11 3). 

This  completes  the  discussion  of  the  algorithms  for  the  computer  program.  A 
few  final  comments  are  in  order  regarding  the  checkout  of  the  program. 


The  program  was  extensively  checked  by  computing  P (or  Q)  for  a large  rante 
of  arguments,  (a,  x),  such  that  both  (hj)  and  (f2)  could  be  used  for  the  same 
arguments.  This  overlapping  procedure  was  also  used  with  tfo)  and  (£j)  , tfo) 
and  @ , © and  @ , © and  © , © and  (fs)  , © and  @ . 
The  subprogram  (ff>)  was  validated  by  comparing  its  results  with  those  from  a double 
precision  version  of  HQ)  . A final  check  that  was  made  compared  the  computation 
of  Q(x  + 1 , x)  from  (w)  , (x  > 100),  with  an  asymptotic  expansion  for  this  quantity 
given  by 


Q(x  + l,x) 


23  23 

+ + 

180  x 2016  x2 


(x*~).  (121) 


Values  of  x as  large  as  106  were  tested  and  the  computed  values  of  Q by  (17)  and 
(121)  met  the  desired  relative  accuracy  for  c = e, , c 1 e2  and  e = eJ(  i.e.,  correct  to 
within  one  unit  in  the  1 2"‘ , 6,h  and  3rd  significant  digit,  respectively.  For  complete- 
ness, the  method  by  which  (121)  was  obtained  and  its  leading  sixteen  terms  are  given 
in  Appendix  C. 

The  program,  as  noted  earlier,  senses  if  for  a given  (a,  x)  P or  Q is  greater  than  the 
specified  e(e,  or  e2  or  e3),  i.e.,  if  (a,  x < is  in  the  zone  of  computation.  If  the  user 
desires  P or  Q computed  for  values  smaller  than  e,  i.e.,  when  (a,  x)  is  outside  the  zone 
of  computation,  this  can  easily  be  accomplished  by  changing  one  location  (ACO)  in 
the  program,  (see  page  H-l). 

A FORTRAN  listing  for  the  entire  program  is  given  in  Appendix  E. 
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COMPUTATION  OFerl(x)ANL)  l/!'(u),  a < 1 


The  formulations  for  the  computation  of  erf(x)  or  erfc(x),  depending  on  x,  are 
taken  from  ( 2 ] . They  are  based  on  rational  Chebyshev  approximations. 

If  ix|  < 0.5 


erftx)  = x Y,  Pk 

k o 

po  = 2.42067  05523  05318  (2> 

P|  - 2.19792  61618  29415  (I  > 
p,  = 6.99638  34886  19136  (0) 
p(  = 3.56098  43701  81539  i 2) 

The  maximum  relative  error  < 10  14,65 
It  1 < x < 4.0 

A 

I erftx t = erfc(x)  = e ' 

Po  = 3.00459  26102  01616  005  C) 
p,  = 4.51918  9537118729  422  (2) 
p2  = 3.39320  81673  43436  870  (2) 
p3  = 1.52989  28504  69404  039  (2) 

p4  = 4.31622  27222  05673  530  (1  ) 
p<  = 7.21175  82508  83093  659(0) 

p6  = 5.64195  51747  89739  71 1 ( 1) 
p7  = 1.36864  85738  27167  067  ( 7) 

The  maximum  relative  error  < 10  1613 


q0  = 2.15058  87586  98012  (2) 
q,  = 9.1  1649  05404  51490(1 ) 
q,  = 1.50827  97630  40779  (1 ) 
q,  =1.0  (0) 


q0  = 3.00459  26095  69832  933  (2) 
q,  - 7.90950  92532  78980  272  (2) 
q:  = 9.3 i 354  09485  06096  2 1 1 (2) 
q3  = 6.38980  26446  5631  1 665  (2) 
q4  = 2.77585  44474  39876  434  (2) 
q5  = 7.70001  52935  22947  295  (1) 
q6  = 1.27827  27319  62942  351  (1  ) 
q7  = 1 .00000  00000  00000  000  (0) 


If  x > 4.0 


erfc(x)  = 


s/n  X 1 k = o 


-2k 


E -k 


,-2k 


(A-3) 


k - o 


p0  = 2.99610  70770  35421  74  ( 3) 
P|  = -4.9473091062  32507  34  ( 2) 
p2  = 2.26956  59353  96869  30  ( 1) 
p3  = 2.78661  30860  96477  88  ( I ) 
P4  = -2.23192  45973  41846  86  (-2) 


q0  = 1.06209  23052  84679  18  (-2) 
q,  = 1.91308  92610  78298  41  (-1) 
q2  = 1.05167  51070  67932  07  (0) 
q3  = 1.98733  20181  71352  56  (0) 
q4  = 1 .00000  00000  00000  00  (0) 


The  maximum  relative  error  < 10'15  61. 

The  reciprocal  of  the  complete  gamma  function  I'(a),  when  a < 1,  is  computed 
from  (25),  (see  also  (1 19)).  The  coefficients  Ck  which  occur  are  given  in  [13].  The 
program  uses  at  most  the  first  seventeen  with  a maximum  relative  error  in  1 /I'(a)  of 
less  than  5 x 10  M.  We  have 


C2  = 
= 


^6  = 

C.- 


1.0 

0.57721  56649  01533 
0.65587  80715  20254 
0.42002  63503  40952  ( 1) 
0.16653  86113  82291 
0.42197  73455  55443  ( I) 
0.96219  71527  87697  ( 2) 
0.72189  43246  66310  ( 2) 


C,  = -0.11651  67591  85907 (-2) 
CJ0  = 0.21524  16741  14951  (-3) 
C, , = 0.12805  02823  881 16  (-3) 
C,2  = 0.20134  85478  07882  ( 4) 
C,3  = -0.12504  93482  14267  (-5) 
C|4  = 0.11330  27231  98170  (-5) 
C1S  = 0.20563  38416  97761  (-6) 
C,  = 0.61 160  95104  48142  (-8) 

1 o 


ci7  = 0.50020  07644  46922  (-8) 


It  was  stated  on  page  38  that  H was  obtained  during  the  computation  for  l/r(a). 
We  have  for 


0 < a < 1/2,  H=  £Ck 


,k  - 1 


A-2 


* . 


which  is  terminated  when  k = 17  or  before  if  k = n < 17  with 


C a"  1 

n 

l 

< - 6 

2 

Then 


1 /I’ ( a ) = a(l  + H|.  0 < a < 1/2  . 


for 


- < a<  1.  H = £ (k(a  l)k  H = — ( I a 
^ (1 


where  H is  terminated  wnen  k - 1 7 or  before  it  k - n < 17  with 

1 


C (a  D" 


< - t 


k - 1 


Then 


/1(a)  = 1 + H - < a < 1 . 


+ fi)  . 


COEFFICIENTS  IN  TAKENAGA’S  METHOD 


In  this  Appendix,  we  list  the  coefficients  associated  with  the  modified  Takenaga 
procedure  as  used  in  (H)  (see  (17)  or  (1 15)).  Also  the  coefficients  for  (18)  or  (80) 
are  listed  as  well  as  those  needed  in  (63).  For  (18)  and  (63)  additional  coefficients 
beyond  those  actually  stored  in  the  program  are  given.  All  of  the  coefficients  were 
obtained  initially  in  rational  form  by  A.  Morris,  (5). 

The  factor  C,  which  appears  in  (17)  is  found  from  e^n  c,  where  the  expansion 
of  Cn  C in  terms  of  R = b is  given  by  (63).  The  expansion  of  fn  C to  b'1 1 is 

fnCs  1/36 b (b  = a - 1)  (B-l) 

-1/81  b2 
+ 1 /360b  3 
+ 2/1215  b4 
69 1/306 1 80  b5 
+ 16/ 15309  b6 
-373/3674160  b7 
+ 80/1 77147  b8 
44069/38972340  b9 
+ 1 792/9743085  b10 
+ 1 1 4953857/63836692920bl  1 . 

From  the  expression  (B-l ),  we  obtain  C in  terms  of  inverse  powers  of  b (see  (80))  as  far 
as  th.-  eleventh  power,  i.e.. 


C = e 


fn  l — 


( >.  i) 


(B-2) 


(68). 


I lie 


j = 0 


j = : 


.1  = 3 

4 


+ 1/36  h 
31/2592  b2 
+ 34 13/ 1399680  b3 
+ 361733/201553920  b4 
1 1388828 1/5079 1587840  b5 
+ 7565202533/78364 1 6409600  b6 
81332300683/1974776935219200  b7 
+ 24588690647475 7/568735 75^343 1 29600  bB 
-1 1336175048324863533/10134871 195854569472000b9 
+ 3090206 1 509879955337/2043 19003308428 120555520b'0 
+ 43461 561 64844627985  I 30999/2390532338'7086090l0499584000b" 

We  proceed  to  give  the  coefficients  NJk  )t  M,k,,  ( for  (17)  or  (1  15).  Also  see 
We  have  from  (115) 


T(a.  x)  = Z(s)  (' 


12  / 

2k 

\ 

h„  * Z « ;k 

L 

N b 

2k. i )/ 

k 1 \ 

) k • 

1 / 

X>  t 

k i \ i k . i / 


a . h are  defined  by  ( 109)  and  (89).  Then  listing  the  N first 


i j 


N 


:k.) 


2k  = 0 

1 00000000000000 
2k  = 2 

0.833333333333333  01 
2k  - 4 

0.555555555555556!)  01 
0.3472222222222221)  02 


- k . ) 

2k  = 6 

j = 4 0.4 1 66666666666  >71)  01 

5 0.6851851851851851)  02 

6 0.9645061  /2839506I)  -04 

2k  = 8 

j = 5 0 333333333333333D  01 

6 0.8 1 900352733686 ID  02 

7 0.3780864197530861)  03 

8 0.2009387860082301)  05 


B-2 


N 

L’k., 

(continued) 

2k  = 10 

2k  - 18 

0.27777777777"’778i) 

01 

j = io 

0. 1 66666666666667D 

01 

0.86955  1 > 243 1 34291) 

02 

1 1 

0.8392801 741 73  22  2D 

02 

0.66!  2838036449  1 SI) 

03 

12 

0. 1 40822229954083D 

02 

0.1 30744 1700960221) 

04 

13 

0.1048597493867351) 

03 

0.3  3489  79  76080384 1) 

07 

14 

0.3  77869445  383984D 

05 

2k  - 12 

0.2380952380952381) 

15 

0.65  7099645  75  785  3D 

07 

01 

16 

0.5048350  73886 1 25D 

09 

0.88 1 85920':  7C934I) 

02 

17 

0. 1 34 1 2654 1 2 1 8996D 

11 

0.90642  292  58  232  "’tM) 

03 

18 

0.534079308498 167D 

15 

0.3151  1 182I4~"502I) 

04 

2k  = 20 

0.32590736396890^1) 

06 

j=  H 

0.1 5 1 5 1 5 1 5 1 5 1 5 1 521) 

01 

0.465 136078722756D 

09 

12 

0.8 16854779541 8281) 

02 

2k  - 14 

13 

0. 1 5 1 645 284022835 D 

02 

0.2083333333333331) 

01 

14 

0.130144560538248D 

03 

0.87  54499958  203601) 

02 

15 

0.571479941  762549D 

05 

0.1  109244092792371) 

02 

16 

0.131 58765492  7 23  2D 

06 

0.5433841632217591) 

04 

17 

0. 1 53805654474373D 

08 

0.1036016990958481) 

05 

18 

0.821 35074 1480055D 

I 1 

0.63258506  ■’062948!) 

08 

19 

0 15509663 1 I87868D 

13 

0.55373342‘70509OOD 

1 1 

20 

0.44506609041 5 139D 

17 

2k  = 1 6 

0. 1 85 1 85 1 8 5 1 8 5 1 85  D 

01 

0. 85966425 54975 89 D 

02 

0.1274420646389341) 

02 

0.79284621 92689421) 

04 

0.2198332488955040 

05 

0.2564060419025480 

07 

0.1002626658580160 

09 

0.5768056531780200 

13 

Nu  . (continued) 


I 


2k  = 22 

2k  = 24 

12 

0. 1 38888888888889D  01 

j = 13 

0.1 28205 1 28205 1 28D  -01 

13 

0.793794836275 189D  -02 

14 

0.7708828747699 10D  02 

14 

0.1 60397649389 170D  02 

15 

-0. 1 67472688827663D  -02 

15 

0.1 545978095 11729D  -03 

16 

0.1 77919959495527D  -03 

16 

0.793808101 128600D  05 

17 

0. 1 038384 1 7905802D -04 

17 

0.22601 7561 582738D  06 

18 

0. 3498760792 28262D  06 

18 

0.355  250664 18604  ID  08 

19 

0.68844702707 1 3*"1D  08 

19 

0.2935 1251 3022 109D -10 

20 

0.7749450300607 15D  10 

20 

0.1 135 1 5291 507799D- 12 

21 

0.4703273760287 15D  12 

21 

0.1 57850106733903D -15 

22 

0.1 36082899 158372D  14 

22 

0.3371 71 28061 7530D -19 

23 

0. 1 434 1 0 1 84689323D  -17 

24 

0.2341 4672265 1062D  21 

1 

2k  + 1 = 3 

ON 

11 

+ 
r i 

2 

0.666666666666667D  -01 

')=  -■ 

0.303030303030303D  -01 

6 

0.85097001 7636684D -02 

2k  + 1 = 5 

7 

0.5  2336860670 1940D-03 

3 

0.4761 90476 190476D  -01 

8 

-0.6430041 15226337D- 05 

4 

-0.55  55  5 55  5 55555  56D -02 

2k  + 1 = 1 1 

2k  + 1 = 7 

j = 6 

0. 2564 1 02564 10256D -01 

4 

0.370370370370370D-01 

7 

-0.87892 149003260  ID -02 

5 

-0.767 195767 195767D-02 

8 

0.789241 622574956D-03 

6 

0.23 1 48 1 48 1 48 1 48 1 D -03 

9 

-0.2 1 568 1 95 1 793063D  -04 

10 

0. 1339591 90672 154D  06 

B -4 


> . 


M2ltt|  . (continued) 


2k  + 1 = 13 

2k  + 1 = 19 

n wnimnimn _.m 

j = io 

0.1 537301 587301 59D-01 

0.880261 7 13595047D -02 

11 

-0.82821 56 153673 18D-02 

0. 1 0 1 29 1 320868040D  -02 

12 

0. 1 4652069306 1 439D-02 

0.425382944673068D-  04 

13 

-0.1 175792689589 14D-03 

0.6243774087 13828D  06 

14 

0.47065769907995 4D -05 

0.22326531 7786923D-08 

15 

0.95201 2948298324D-07 

2k  + 1 = 1 5 

16 

0.92092730290985 ID  09 

0.19607843 1 372549D  01 

17 

-0.364022597059366D-1 1 

0.868360035026702D  02 

18 

0.3845371021 18680D-14 

0.1 1961 470270395  7D -02 

2k  + 1 = 21 

0.66655778057 1596D  04 

j=  11 

0.1 4492753623 1884D-01 

0.1 561 75357088423D  05 

12 

-0.8053446095 24 176D  -02 

0.1 3799922975591  7D  07 

13 

0.1 56253676 162327D -02 

0.3 100907 19 148504D  10 

14 

0.1 4249827604399 ID  03 

2k  + 1 = 17 

15 

16 

0.679474426392697D  -05 
-0.1 75  14133195771  5D-06 

0. 17543859649 1 2 28  D —0 1 

0.84985 142291 398 1 D - 02 

17 

0.240079 1 0704 1 6 1 2D  - 08 

0.1  344860333506720  02 

18 

-0.1 6288 163096471 3D  -10 

0.9206186407701 27D  04 

19 

0.4622502 19327808D  13 

20 

0.3560528723321 1 ID  16 

0. 2939823607353 74D  05 
0.4 '’70763 58950073 L)  07 
0.2454 1 4654868959D  09 
0.369 1556 1 8033933D  12 


I 


r - wrr-- ' - 


MMiip.l  uj.  I IW  I p pi  |W  JP 


(continued) 


2k  + 1 = 23 

2k  + 1 = 25 

j=  12 

0.133333333333333D 

01 

J=  13 

0.1 234567901 23457D 

01 

13 

0.78228731 5947945D 

02 

14 

0.7596263771 72480D 

02 

14 

0 1 641 235305541 89D 

02 

15 

0.1704820355961030 

02 

15 

0.16641 2253953536D 

03 

16 

0.189106579970833D 

03 

16 

0.9136909627549970 

05 

17 

0.1  16720392056449D 

04 

17 

0.28427 1767098726D 

06 

18 

0.4227344701  27948D 

06 

18 

0.5034342786294050 

08 

19 

0.9 1 3761 861 504362  D 

08 

19 

0.49 1 1 82369 1 3403 1 L> 

10 

20 

0.1 165340093253870 

09 

20 

0.243  P07 1607648 10 

12 

21 

0.8401 5 1056535048D 

12 

21 

0.5 1 25 1 267737405 1 0 

15 

0.3131852910855800 

14 

0.2967107269434260 

18 

23 

0 5038430'725  144080 

17 

24 

0.2247808537450200 

20 
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ASYMPTOTIC  EXPANSION  FOR  Q(x  + 1,  x)  (x  * ~) 


In  this  appendix  we  show  how  the  asymptotic  expansion  of  (121)  was  obtained. 
The  form  of  the  expansion  is  given  by 


1 1 [2  M " 1 

Q(x  + 1,  x)  a - + - — Y\  e /xk, 

2 3 Vrrx  k“0 


(C-l) 


where 


e0  = l,  e,  = -23/180,  e2  = 23/2016, 


as  shown  on  page  43.  The  ek  for  k = 0,  1,  • • • , M - 15  are  given  on  page  C-8  of  this 
appendix. 


We  start  with 


Q(x  + 1 , x)  = r(x  + 1 , x)/r(x  + 1 ) (See  (2),  (4),  (5)),  (C-2) 


where  the  asymptotic  expansions  for  the  two  functions  on  the  right  hand  side  are  given 
by 


. -i  m - 1 

1/P(x  + 1)  = l/x  T(x)  = e-  fn  * r(x)  ex  (n/2tT  xx  + ,/2)  ^ fk /jck , (C-3) 

k = 0 


with 


1 


B, 


fnxP(x)a(x  + -Wnx-x  + -fn2n  + ^ , (x  + °°),  (C-4) 

\ 2 / 2 , 2m(2m-l)x2ra-1 


(B2m , the  2m,h  Bernoulli  number),  (See  [ 1 ] , p.  257,  8 1 0), 


and 


2 M 


r(x  + 1 , x)  =£  e~x  xx  + 1/2  £ bk  r(k/2)/x(k  0/2 , (x -*•  °°). 


(C-5) 


k = 1 


C-l 


The  coefficients  bk  r (^J  and  fk  are  constants.  The  multiplication  of  the  series 
in  (C-3)  and  (C-5)  to  yield  (C-l)  was  carried  out  by  A.  Morris  using  his  “Flap’  routine, 
[5J.  In  addition,  he  also  used  “Flap”  before  hand  to  determine  the  fk  and  bk  r 
in  rational  form.  For  completeness  they  are  included  on  pages  C-6  and  C-7,  respectively. 

The  bk  were  computed  independently  by  M.  Saizan  in  decimal  form  and  are  given 
below. 


b,  = 
b2  = 
b3  " 
b4  = 
b5  = 
b6  * ' 

b7  = 
b8  = ■ 

b9 

b.o  = 
b.,= 
b,2  = 
b,3  = 
bM  = 
b.S  = 
b.6  " 
b.7  = 
b.8  = 
b I 9 = 
b20  = 
b3.= 
b32  = 
b23  = 
b24  ~ 
b25  = 
b26  = 
b27  = 
b28  = 
b29  = 
b30  = 


707106761186547 524400a 443621 D*QQ 
6666666666666666666666 666667 0+00 
117 8 511 30 19 77 57 92 073347406640+00 
2 962962  962962  962962  962  962  9620*01 
32736426054932755759298350060-02 
1410 934 7442 68 0776 0 1410 9 34 7450-U2 
1 01 1191 7961 412562 33453 3o2 3690  Oc 
3 1354 10  542  81 7 950  2 2535  76  82  762  D-0  3 
24725527389373814204443595020-04 
29664995371442559371228781380-04 
18773314722352835450708456900-04 
56531048757843453773952173490-05 
30356279857570629915558755250-06 
65675582619137471472473326040-06 
3966166  25158  6 64330 109 10  363480- 06 
11709055465263091499753584840-06 
4618562455125  90  916288  4910  668  0-00 
14926776659329080172708068980-07 
8819467  30  7943  99  70  6342  56  29164  0-08 
2 5741 6667141 8 45 790 161 43ol 5190-08 
7 9680 30 9541 18 5471 3047 7 3 765 7b0-10 
3 452685 580 6 98 60 90 2 791 45b36 11 0-09 
20163942557072502446885174830-09 
58439462516833163681724231240-10 
14884457788493658882335740850-11 
8 09  05  37  2853  65  5 317  8537  829996  0 0-11 
4691 7438130337419702825*47300-11 
13535257572473374162248894010-11 
29316941779385508468669639750-13 
19147882067656010283127746780-12 


We  note  the  first  three  terms  of  (C-5)  are  given  in  [1,  p.  263] . 


•The  symbols  bk  and  aR  are  used  in  this  appendix  to  retain  the  notation  in  (6). 
They  should  not  be  confused  with  those  in  (109). 


C-2 


\ •* 

» . 


The  derivation  below  for  (C-5)  is  given  for  completeness  and  is  based  on  the  mate- 
rial given  in  (6,  p.  85  -87) . 


We  have 


l'(x  + 1,  x)=  f tx  e * 
Jx 


Letting  t = Ml  + u),  (C-6)  transforms  to 


r ( x + l , x) 


= e ' xx"f  e x,u  fn(" 


u)l  du. 


The  integral  on  the  right  in  (C-7)  is  treated  in  |6.  p.  87).  Its  asymptotic  expansion 
(x  * °°)  is  obtained  by  a series  reversion.  Indeed,  letting 


, u2  u’ 

V = u i n( ! + u I = — — + 

: 3 


. lu|  < 1. 


we  note  that 


d\  , 1 _ u 

dll  I + U 1 + U 


du  1 + u 


d\  u 


Then,  assuming  an  expansion  for  u of  the  form 


u = \3v  + £ ak  Yk/:.  (a,  = VT  ). 


we  have  b>  differentiating  (C-10)  and  using  (( -CM 


1 + u * k k i _ 

— = E -ak  v’  • ai = 

u k i - 


K" 1 1 


It  follows  from  this  result  that 

OO  w OO  L OO  I L 

•■*£•.*>-£  •.  v!  E ',vl" 

k 1 k = I k = i ^ 

Expressing  the  product  on  the  right  as  a Cauchy  product,  (C-ll)  becomes 

OO  ^ OO  | i 

1 + E ak  yI  = E yI " ' Z t s 3m 

k - 1 i = I j = I - 

which  contains  a Unear  set  of  equations  for  the  ak , k = 0,  1 , • • • with  a0 
example,  equating  coefficients  of  like  powers  of  V,  we  have 

1 „ r~ 

0 = - a,  aQ  =>  a0  = 0,  since  a,  = V2 

1 , rr 

1 = - a:  + a,  a =>  a = \J2 

->  I l 0 1 

1 

a.  - ; ai  a2  +a!  a2  =>  a2 


V 1/2 : 
V° : 
V,/2 : 


In  general,  by  setting  k = n 1,  i = n + 1,  (C-12)  can  be  solved  for  an. 


" i 

a , = V*  - a • a . (a„  = 0) 

n I j * I j 0 


1 I “ 


or 


a = 


n ( n + 1 )a , 


n 1 i 

y - a a 

i-J  T J '•  * 1 I 


) ■ 2 - 


n = 3.  4, 


Thus. 


a,  = 2/3.  a_j  = a4  = 2/135,  a^x/2/1080.  aft  = 4/8505, 


etc. 


1 

(C-ll)  ‘ \ 


(C-12) 


0.  For 


We  get 


(C-13) 


C4 





The  expression  (C-13)  gives  the  coefficients  for  the  inverted  series  in  (C-10).  Therefore, 
following  [6] , the  coefficients  bk  of  asymptotic  series  (C-5)  can  now  be  determined. 
We  have 


! 

s 


du 

dV 


v(k  2)/2 


k I 


The  bk  are  the  desired  coefficients  that  appear  in  (C-5).  Thus,  e.g., 


\ls/2,  b2  = a2  = 2/3,  b3 


b4  = 2a4  = - 4/135. 


The  ak  are  listed  below  for  completeness. 


a = . 1414213562373Q9504880163S724D«-01 

a'  = . 66b66666666666666666666666670*00 

a2  = » 78567420131838613822316040230-01 

a4  - .14414614614614814814814814810-01 

a = . 13094570021973102303719340030-02 

a = .47031158142269253380364491490-03 

a = .28891:94175464463612968067690-03 

a = . 78385263570448755633940  319050-04 

a9  = .54945616420830698232096877820-05 

a*  = . 59329990742  885118742457562760-05 

a ° = .34133299495186973546742648900-05 

a = .94218414596405756289920289150-06 

a = .46  701969  01164  712294  7013469610-0  7 

a = . 93822260884482816389247609770-07 

a'  = . 52862216667621911747880484620-07 

a'5  = . 146  36  319  33157  88643  74  691  9810  50-0  7 

a'6  = . 54336028883834225445704831380-09 

a'7  = . 1656530  73992 545424 1 41  20  0 76650-0  8 

a‘“  = . 92841761136252600667638201730-09 

a!’  = . 257416667141845  79016143615190-09 

a = .7588600 9086843306004546443550-11 

a2'  = . 31388050733623718435632396460-1  0 

a” = . 17533  86  3 09310  65238  66856678110-1  0 

a = .48  6 99  55  2097  36  09  69  734  7 701927  00-11 

a2  = . 11907566230  7949  27 1 0 58  6859268  0- 1 2 

a = .62234902195119475272140  768920-12 

a = . 34  753  65  78  74  32  4014  594685  8868  90-12 

a2(J  = .96680411231952672587492100070-13 

a = . 20  218  58053  750  72 472 1 992388 9480-14 

a = . 12765254711770673522085164520-13 


C-5 


The  rational  coefficients  f for  (C-3)  are  included  in  the  expression  below,  (C-14). 
The  quantity  T = 1/Vx.  The  coefficient  of  1 /x3 , f3 , is  given  by  1 39/5 1 840. 


-1/12*T**2 


1/288*T**4 


139/51840*T**6 


-57 1/2488 320* T* *8 
-16387 9/2090 18880 *T* *10 


5246819/75246796800*T**12 


5 34703 5 3 1/902 96 156 1600 *T** 14 


-448 3 13 1259/86684 309913600*T** 16 


(C-14) 


-432261921612371/514904800886784000*T**18 
6232523202521089/86504006548979712000*T**20 
2 58 3 462 9665 134204969/1 34946 2 502 16408 3507 2000 *T** 22 
- 1579029 1388 54 91 9086429/97 161 300 1558 1401 25 1840000 *T**24 
-74659086996265 1602 20 3 151/1 16 59 3 5601869768 1502 2080000 *T** 26 
15'1513601028097903631961/2798245444487443560529920000*T**28 


The  coefficients  bk  I'( k/2 ) for  (C-5)  are  included  in  the  expression  below,  (C-15). 
The  quantity  T = I ly/x.  The  coefficient  of  l/x3/2,b4  r(2),  is  given  by  - 4/135.  The 
coefficient  of  1/x3,  b7  1(7/2),  is  given  by  ( 139/103680)  \/Tn. 


1/242«*11/2>=PI  ->=11/21 


2/3-1 

1/254/  •>=(  i /2  1 ->P  l - i 11  /2  1 ■>  ! = ■>/ 

-9/  1 35  = I »*3 

1 /5  76»2  = = ( 1 li  I “PI  »»(  1 / 2 I ( **/. 
fa/2o35-'i  ■'  ■>  5 

-13-»/10jSor  = 2 = =ll/2  1fP:->tll/2)  r - - 6 
It/a5j5->l  = *7 

-b71/’s976t^-i>2  = = ll/2l-PI<-;(l//l-I-=6 
-u 9 92/1^6/ 9, 23  f » - 9 

lo3o79/9lc.37/V*  /4-ii//l-Pl4-ti//l  - T - ■>  1 7 
-3391*9/952367075-14411 

5<5b819/l5C.<.9  35V1oOU'>/  = =U/^|iPI*i(l/2l=T-=)2 
636  152/1977  / ni?25*T  = = j 3 

-,397*  35  3 1 /Iff,  5 >2  312’ 2 00-  2 »' s 1 1 /2  ) - PI  - (1  //  I -T  19 

(C-l  5) 

23795C  12225/3^5oc95lJ2'/9373=I  = ->li 

-9983131259/!77258fl9S27?5C-2-  (l/2)-‘r-l>  (l/2)-l»->lb 
-13573r529713h/2?5527Cfib70f--7  7 5-  T = 917 

<u2261?clbl2771/lC298  0 9fcC1773  36of.,0-2ioil/<'i«hl-9(l/2)'. 
-o3195259/3392/676569«0oll9il*5'-T*®19 

6/3/5/3/72521065/17  30  *601309795  55 1 4*  r,i)v24i<l/2|..Pi4,M]/2»»T'S-*«0 
a773995C«t01«olb/ 7G029912/U758  0937-, -' !->-?) 

-/  5c39<-/9t631  39^099  t.o  / 2b9bQi  5*093*  M b 70  )o  * ■ . u 0 ~i  • * l 1 / 2 ) - -H  I * = l 1 / 2 ) - 1 2/ 

930*997  /C/2t59969/210*7973frb3/f"--*b331?54  I - 23 

-1579*27l3p955919j3  65/9/195j2/6*03li62ll>.23f.3()Olib0  0;i':-(l/2l-P|-':(l/^)-7-'25 
-i7*5o5,*  7 9 237n038/P9/*<,il3f55„5jri9/715<-2->-i  -2» 

75>-59*P097f/b3]  b0t2031*;l/^33lo7i2*3?7  353b3rj5«,lb0*tJ'*'^->'-l  1 / 2 i ■ P 1 >-(  4 / * ) - I - -4 
-•,6o3b0  37i.-  35  3b  70  3K2b*9/P/6c97r;/7,7-,9o2P7/8f7*171c73  ! 4-27 

1 3l  1 5 4 3bO  1 T *6*9  79*3b31  9bl  / 5-,5b5  y'.c  Ho9  /5r.P  M?  1 *5  93*0*0'*  * - - 11  / * ) > . - ( 1 / * I T- 

?o5  311733. -P  02  60  667  3693900  6/3  3*0  5 6*2/  3 o5  5 57-153-9^5,  n5j  7 5-14  4/  5 
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The  coefficients  eR  for  (C-l ) are  included  in  rational  form  in  (C-16)  below  and 
then  ek , k ^ 0,  are  given  in  decimal  form  in  (C-l  7).  Again  T = 1 !\f\.  The  coefficient 
of  1/V’/2  from  (C-16)  is  ( 22>l21Q)l\f2n . 

1/2 

2»9/77  7*P*?»M-l/fl»PI»»(-J/>)»I»»7 
-2016j7J/lf<.7S7  7AOC«P'»(-J/il«Pi*»(-l/?l»r»»9 
-*.568j87/i.C75  86»)*C*?»»I-1/<I«PI»»(-1/?IM»»I1 
*7839y$9/l 1 id 70  19S7&(W»«  1-1/2  )"P» •«  I - l/?l  ® 1 • • 1 i 

1 132  1 1 961  792«7/12Cl88  79  790S92Cil00»2*»  I - I /?  I*PI  ••  l-l  /2  I*  S , 

((.  -16) 

-:7<.0(,704t0353/2<.7C5C28i2S3  7M,0»7>'»l-l/2»»PI»»l-|/7)»I*»17 

- 1)U 3 30S782  3 740  1 1 /772  3 5 720  l 334) 76C0J»2*»< -1 /2  1»PI •» i-1  /}  I* 7 ••  19 

62  8 7 3i4d424375257233/‘.4  76>82  3383n‘y7v,n0V8nil0()*2o»(  -1/2  I •PI»»t-l/2l»I»»21 

14145  / 3 on 2 4 90 32 U7 3 33/4048  3 8 7S06-.9,;7S0S7 1800 •>••(-!  /2  l«P|»»  1-1/2  )»!»»23 

-4  32  054y0  7ClloRi91445  1 5»3  / 1 C20  1 93t5 1 6 3RC4  7 1 3 1 44  32  ut)C0°2»»  1-1//)«M*»(-1/2)»T»»2  5 

-93i2/9S3595952  73C,C2  7426859L0  7/753  1 o52504  1 78734#0»3«B«1280000u»2»»l-l/2  l*PI  ••  1-1 /2  1 •! »»27 

23 3460 96240 7729 55 7 3*46030*47/ 1301 1141 3 1 alRtbl 253*4*4 1 2l0000»2»» l-l /2 I *P1»» 1-1/2 )*I»»29 


Decimal  form  lor  e.  , k > 0. 

k * 


c0 

= 

1 .00000 

00000 

00000 

00000 

00000 

00000 F 

00 

c, 

= 

-1.27777 

77777 

77777 

77777 

77777 

77777F- 

-01 

C2 

= 

1.14087 

30158 

73015 

87301 

58730 

15873F 

02 

CJ 

= 

4.99614 

19753 

08641 

97530 

86419 

753081 

03 

C4 

= 

1.63704 

05768 

07166 

31333 

91528 

45 325 F 

03 

C5 

= 

1.68125 

66626 

29637 

16667 

42037 

1 1 240F 

03 

C* 

= 

9.01554 

11107 

14762 

01363 

03263 

87528F 

04 

C7 

= 

1 .40480 

10663 

69107 

331  17 

86737 

7 6069 F 

03 

CH 

= 

1.05651 

27415 

16873 

01128 

361 10 

5 6040 F 

03 

S 

= 

2.15439 

00162 

97250 

29517 

66213 

759 I6F 

03 

cio 

= 

2.10674 

16054 

60398 

61407 

67106 

9I024F 

03 

C|. 

= 

5.24124 

92389 

44466 

72838 

31327 

57I08F 

03 

e.2 

= 

6.35254 

25734 

46898 

69625 

71691 

5 4688 F 

03 

e.3 

= 

1.85264 

1 1158 

31  105 

51190 

42883 

84580F 

02 

(u 

= 

2.69132 

88833 

12297 

39720 

54734 

83920F 

02 
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CRITERION  FOR  TERMINATING  THE  SERIES  IN  (10) 


In  the  main  text  of  this  report  an  expression  lor  P(a,  x)  is  given  by  (10),  i.e.. 


I’la,  x) 


e ' x" 


al  l a ) 


I + 


a + I (a  + 1 )(a  + 2 ) 


+ s + I m 

ii  N ♦ I 


(D-D 


where 


s x"  / 1 ( a + I )(a  + 2)  • ■ • (a  + n)| 


T , = s , + s , + • 

n * I n ♦ I n ♦ 2 


(D-2) 

(D-3) 


The  objective  in  this  Appendix  is  to  prove  the  statement  made  in  Section  5 
dealing  with  (l)0).  Given  an  c > 0,  if  N is  defined  as  the  smallest  positive  integer  for 


which  sN  < — . we  wish  to  show  the  truncation  error.  TN  t , , associated  with  the  series 


in  square  brackets  of  (D-l ) is  less  than  e . provided  x/(a  + N + 1 ) < 2/3.  The  proof 
requires  only  a few  lines. 


From  (D-3) 


1 n . I SN 


l * \ 

X X2 

] + + + ■ • • 

\a  + N + 1 / 

a + N + 2 (a  + N + 2)(a  + N + 3) 

< 


- • - • ( 2/3  )k  = c 


where  we  have  used  the  inequalities  x/(a  + N + j)  < x/( a + N + I ) < 2/3,  j = 2,  3.  • 


!)-l 


1 

1 
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FORTRAN  LISTING  OF  THE  PROGRAM 


I 


4 


I 

i 

« 

? 


K 

\ 

i 

i 


I 

i 


Calling  sequence  to  the  incomplete  Gamma  Function  Ratio  Subroutine. 
CALL  PAX(a,  x,  P,  Q,  10P)  where 


IOP  = 0 for  12  significant  digits  of  accuracy  with  an  error  no  greater  than  1 unit 
in  the  12,h  significant  digit  of  P and  Q. 

IOP  = 1 for  6 significant  digits  of  accuracy  with  an  error  no  greater  than  1 unit 
in  the  6Ih  significant  digit  of  P and  Q. 

IOP  = 3 for  3 significant  digits  of  accuracy  with  an  error  no  greater  than  1 unit 
in  the  3rd  significant  digit  of  P and  Q. 

The  routine  is  designed  in  most  cases  to  set  P = 0.  0 = I if  P < c and  P = 1,0  = 0 
if  0 < e.  The  quantity  e is  ordinarily  set  to  take  the  specified  value  e = 5 x 10  ' 
5x10  7 , or  5 x 10  4 when  IOP  is  set  to  0.  1 , 3,  respectively.  However,  c may  be  set 
internally  to  zero  or  any  other  positive  value  less  than  the  value  of  c specified  by  IOP. 
In  the  program  c is  stored  in  ACO. 

Restrictions  x > 0 and  a > 0. 

I.rror  Returns:  If  either  a < 0 and/or  x < 0 then  P is  set  equal  to  2. 
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SUBROUTINE  PAX  ( A , X , ANS , Q ANS  , I OP3 > 

OIMENSI ON  VA(2C)«V3(20I  tZA(30)»ZB(30)  ,S A 120 ) ,S8  (20  I 
DATA 

P3/2. 4266  79552  33532  E*2 


1 

2 Pl/2.1979  26161  82942  E*1 

3 P2/6.9963  €348*  61914  E»0 

4 P3/-3.5609  84370  18154  £-2 

5 Q3/2.15C5  88758  69861  E + 2 

6 Ql/ 9.11 64  9054C  45149  £»1 

7 Q 2/1. 50  82  79  763  C4078  E ♦ \ 

8 Q3/1. 

9 CO/3. 0045  92610  20162  E+2 

OATA 

1 Cl/4.5191  89537  11873  E*2 

2 C2/3.3932  08167  34344E*2 

3 C3/1.5298  92853  46940  E*2 

4 C4/4.3162  22722 

5 C5/7.2117  582  3'' 

6 C6/5.6419  55174  7897*  E-l 

7 C7/-1.3686  4857  3 82717E-7 


20567  E*1 
8 8 3 'J  9 E»0 


/» 

/. 

/. 

/. 

/. 

/, 

/. 

/. 

/ 

/, 

/. 

/. 

/» 

/» 

/. 

/. 


9 ALPH3/-2.9961  C7C77  C3542  E-3/, 

9 ALPH1/-4. 947 3 09106  23251  £-2  / 

OATA 

1 ALPH2/-2 .2695  65935  39687  £-1  /, 

2 ALPH3/-2. 7666  130*6  C9648  E-l  /, 

3 ALPH4/-2.2319  24597  34185  t- 2 /, 

4 BET  0/1. Cb20  92335  28468  £-2  /, 

5 OCT  1/1. 91 jC  89261  07*30  E-l  /, 

6 BET2/1.C516  751"7  36793  /, 

7 BET 3/1.  9873  32018  17135  /, 

8 BET 4/1.  /, 

9 03/3.3045  92609  56983  E*2  / 

OATA 

1 01/7.9095  09253  27898  E*2  /, 

2 02/9.3135  43  948  5361'  Et2  /, 

3 03/E.3698  02644  65631  E *2  /, 

4 04/2.7758  54447  43988  E»2  /, 

5 05/7. 7CCJ  15293  52295  E»1  /. 

6 06/1.2782  72731  96294  £*1  /, 

7 07/1.  / 

CST2  = SP< OP 3413. /OP  13996  80  . ) 
CST3=SP(OP361733./OP2; 155  3523.  » 

CST4=SP(OP 111383  8 2ol . /0P5  0 7 91  58  784  0.  ) 
CST5=SP(OP75652  0 2533. /0P78 364  16439  €00.  I 
OATA 

1 CST5/. 96539  05736  46935E-3  /, 

2 CST4/. 22422  666CC  5CC12E-2  /, 

3 CST3/. 17947  20737  75593E-2  /, 

4 CST2/. 24334  14494  74166E-2  /, 

5 GANPT5/1.7724  5385C  90552  /, 

6 RT2PI/2.5C66  28274  63130  / 

7.RTPI/1. 7724  53350  90552  / 

OATA  ACC1/5.F-1 3/. ACC  3/ 5. £-7/ , ALPHA  3/4.892/. 

1 ALPHA  1/7. 1316/,  ACO/5.E- 1 3/» 

1 XOl/31.CC0/,X03/l7.C0/*CEP/2.3O25  65092  99405/ 


i.-: 


2 ,UXl/l.i366/,UX3/2.9245/ 

3 ,RT2/. 14142  13562  3731GEU / , ABAR/130 ./ 

4 ,ALPHA5/3.29053/,X35/9.7/ , UX5/3 . 8 87/ , ACC5/5. E-4/ 
ASA V=A 

I NO  = 0 
ANS=0. 

IF  ( IOP3.NE.O  ) GO  TO  1105 
ALPHA=ALPHA1 
X0=X31 
ACC=ACC1 
UX=UX1 
GO  TO  1111 
1105  CONTINUE 

IF  ( IOP3.NE.1  > GO  TO  1147 
ALPHA=ALPH A3 
ACC- ACC3 
X0=X03 
U.'=UX3 
GO  TO  1111 
1137  CONTINUE 

ALPHA=ALPHA5 
ACC= ACC5 
XC=X05 
UX=UX5 

1111  CONTINUE 

IF  ( A.LE.C.O  .OR. X. IT. 0.0  ) GO  TO  1131 
IF  { X.GT.C.O  » GO  TO  1151 
1121  CONTINUE 

GO  TO  1171 
1131  CONTINUE 
ANS=2. 

RETURN 

1151  CONTINUE 

IF  i A. IT. 15.  J GO  TO  1331 
RTA=SQRT (A) 

IF  C ACO.LT.ACCl  I GO  TO  1155 
IF  ( A.GE.X  I GO  TO  1161 
FP=1.-1 • / I 9. * A) ♦AlPHA/(3  ,*RT  A ) 

FP=-FP*FP*FP* A»X 
IF  ( FP.GE.O.G  > GO  TO  1191 
1155  CONTINUE 

IF  ( A.LT.ABAR  » GO  TO  1331 

IF  ( <.75*A).GT.X.OR.(l.25*Al .IT. X ) GO  TO  1331 
GO  TO  5 U 11 
1161  CONTINUE 

IF  ( X.LE.UX  I GO  TO  1171 
FN=1.-1./(9.*A)-ALPHA/(3.»RT A) 

FN=-FM*FM*FM*  A»X 
IF  < FH.Lt. 0. C ) GO  TO  1121 
GO  TO  1155 
1171  CONTINUE 
ANS=Q. 

QANS-1. 

RETURN 

1191  CONTINUE 
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1311  CONTINUE 
ANS-1. 

Q ANS-C . 

RETURN 

1331  CONTINUE 

IF  ( A.LE.X  t GO  TO  1441 
1341  CONTINUE 

IF  ( A.GE.3C.  I GO  TO  1343 
ALGX=ALOG(X) 

T1=-X»A*ALGX 

CALL  GAMC03  ( A, H, ACC. CANS  ) 

R=EXP  <Tl ) /CANS 
GO  TO  1^49 
1343  CONTINUE 

AISG=1./ <A*A ) 

T1=A-X»A*AlOG<X/AI4 ( ( (AISQ-l. 3333  33333  33333)*AISQ*4.b66  66666 
1 666667  ) VISQ-14C  • I / (IfefiC  . * A 1 

R=RTA*EXP<T1 ) *. 39894  226C4  01433 

1349  CONTINUE 

IF  ( A.GT.X1  GO  TO  1351 

1350  CONTINUE 

IF  t A.LT.l.  ) GO  TO  1459 
IF  C A.GE.X  ) GO  TO  3311 
GO  TO  3005 

1351  CONTINUE 

IF  ( R.LE. (A* U.-X/ <A4l. >1*AC0> » GO  TO  1171 
GO  TO  1 35G 
1441  CONTINUE 

IF  I X.GE.XO  > GO  TO  1341 
TW0As2.*A 
J = INT (T  HOA) 

Tl=ABS(FLOAT (JI-TWOA) 

IF  ( Tl.GT.O.JGO  TO  1341 
I = J/2 

AF= A-FL  OAT ( 1 1 

IF  ( AF.GT.C.  » GO  TO  6011 
GO  TO  6071 
1459  CONTINUE 

IF  < X.LT.1.5)  CO  TO  1461 
Tl=  (2.*X-A )/  (1. ♦*> 

IF  ( R.LE.  (ACO'X’Tl)  ) GO  TO  1311 
IF  t R.GT. (. 101*X*T1 I ) GO  TO  3011 
GO  TO  7011 
1461  CONTINUE 
INO=l 

GO  TO  3C11 

THIS  IS  PART  F 

2011  CONTINUE 

ACC7=.5*ACC 
T 7 = A*  AL  G X 
CEE  =2 . 

T=T7 /CEE 
SUM=T 


I -4 


2031  CONTINUE 

CEE=CEE*1 • 

T= ( T *T7  > /CEE 
SUH=SUM*T 

IF  ( ABS  (T) • GT.ACC7  ) GO  TO  2031 

T2=T7* (1  .♦SUM  I 

TK=-X 

CE£=1. 

A JK=TK/  ( A4-CE E > 

SUM= A JK 

2071  CONTINUE 

C£E=CEE ♦ 1 . 

TK=-(X*TK»/CEE 
A JK=TK/ <A*CEE» 

SUM=SUM»AJK 

IF  < ABS (AJK) .G£«-ACC*SUH)  GO  TO  2071 

T3=-A*SUM 

T 1 = H 

T4=Tl*T2 

T5=T1*T2 

QANS  = T3-T4*T3MT4*T5)-T5 

ANS=1.-QANS 

RETURN 

3005  CONTINUE 

IF  ( X.LE.CEP  I GO  TO  3011 
IF  ( X.LT.X3  J GO  TO  7011 
IF  ( A.LE.2.  I GO  TO  3006 

IF  C (RMX*l.n  .LE.|ACO*X*(2.»X-A>)  > GO  TO  1311 
GO  TO  9C11 
3C0  6 CONTINUE 

IF  < R.LE.  (ACOMX«l.  -AM)  GO  TO  1311 
GO  TO  9011 
C 
C 

C THIS  IS  PART  A 

C 

3C11  CONTINUE 

HAVACC= ACC*. 5 
ROVA=R/A 
SNP1 =1 . 

SUM=SNP1 
CEE  = A 

3031  CONTINUE 

CEE-CEE*!. 

Tl=  X/CEE 
SNPl=SNPl»Tl 
SUM  = SUM  *SNP1 

IF  ( Tl.GT.. 66666  66666  66667  ) GO  TO  3031 
3041  CONTINUE 

IF  ( SNP1.LE.HAVACC  » GO  TO  3051 
CEE=CEE»1. 

SNP1 = ISNF1*X ) /CEE 
SUH  = SUM*  SNP1 
GO  TO  3 C 41 
3051  CONTINUE 
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ANS=ROV  A*SUN 
QANS  = l.-At,S 

IF  ( IN0.E0.1.AN0.ANS.CT.C.9  » GO  TO  2011 
RETURN 

THIS  IS  PART  El 

5011  CONTINUE 
ASA V= A 
A=A"1. 

ALPHY=3 • *SQET ( A ♦ 2 . / 3 . I 
AlPlNV=l./ALPHY 
A2THRU=A*. 66666  66666  66667 
DELTA=(X-t)'. 66666  66666  66667 
SZ=QElT A/A2THRO 

Z=E  XP ( . 3 3333  33333  33333*  ALOG  C X/A2THRDH 
S0= t ALPHY’SZ ) /<  <Z*i.)*Z*l.) 

XSA V=  X 

X= (S3*S0 )/2. 

RTX= ABS (SC  J/RT2 
T5=EXPC-XI 

IF  ( X.GT.0.25  » GO  TO  5361 
Tl=  < <P3*X*P2)  *X*PD  *X*PC 
T 3= I <Q3*X»02»*X»C1)*X*QO 
DEL  = RTX*  < T 1 / T 3 1 
CERF=l.-OEL 
GO  TO  5361 
5361  CONTINUE 

IF  < X.GT.16.)  GO  TO  5371 

Tl-CQ+RTx*  (C1*RTX*(C2*RTX* <C 3 *RT X* ( C4 »RTX * <C5 *RTX* (C6*C7*RTX>1)>  J 

1 > 

T3  = 03*RTX*<D1*RTX*<02  4RTX*<03*RTX*(04*RTXM05*RTX*<06*RTX*D7>  )»l  ) 

1 > 

CERF=T1/T3 
DEL=1.-CERF 
GO  TO  5361 
5371  CONTINUE 
T=1 ./X 

Tl= ALPHC*T*( ALPH1»T* ( ALPH2+T • ( ALPH3 *T *ALPH4» ) ) 
T3=BETC+T*(QET1*T*(0ET2»T*(9ET3*T*0ET4))» 

CERF=  :i  . /~TX  ) Ml  ./RT  PI  ♦ T 1/  ( T 3 *X  I ) 

OEL=l.-CERF 
5381  CONTINUE 
ZA0=*1. 

IF  ( X.LE.0.25  I 230=1 . 5*CE P F *RT2P I ) /T5 
IF  ( X.GT.C.25  ) Z3t=.5*CERF*RT2PI 
IF  ( SC .GC.C.  ) GO  TO  5389 
Z AC  =-l • 

5389  CONTINUE 
X-XSAv 
Zb«l)=Zi 0 
Z A C 1 1 = Z A c 

SCSO=Sl *s 
T 7 = -l . / (A.)S(Sv)  » 

T OL  = ACC  * ZUO 
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00  5391  L=2 » 5 

K=2*<1-1) 

T7=T7*S0SQ 

ZB(LI=-T7 

♦(FLOAT (K'-1.)*Z8(L“1I 
Z A (L ) =-T  7 *S0 ♦FLOAT ( K ) *2A ( L“1 ) 

CONTINUE 
VB (1 )=Z8 (1 ) 

VB(2I  = -.83333  33333  333 33E-1 • ZB (3 1 
VA(1>=. 66666  66666  66667E-1*Z A (3 > 

SB (1 1 =V  8(1 ) 

ALPHSQ* ALPINV*ALPINV 
ALPHCB= ALPHSQ* A LPINV 
S8(2I=V8(2>*AIPHSQ 
SA(1I=VA(U*ALPHC9 
T15=SB(1)*SB(2> 

T 17=SA ( 1 ) 

VB(3)-  -(Z6(4l-.C625*ZB(5))*  .55555  55555  E5556E-1 

VA(2I=. 47619  04761  90 476E-1* Z A (4 )-. 5555 5 55555  55556E-2* Z A (5 > 

T23= ALPHSQ 

15=3 

MF=0 

CONTINUE 

17*15-1 

T23=T23* ALPHSQ 
SB(I5)=va(I5) *T23 
SA(I7)=VA(I7) *T23*ALPINV 
T15=T15*SB(I5I 
Tl7=Tl7*SA(I7) 

IF  ( ABS (S  B ( I 5 ) ).LT.T0L.AN0.AGS(SA(I7I)  .LT.TOL  I GO  TO  5431 
13=2*15 

1 5=15*1 
T7=T7*SQSQ 
T25=4.*FL0AT ( 15 ) -7 • 

ZB(I3)=-T7+T25*ZB(I3-1J 
ZA(I3)=-T7*S0 ♦(T25*1.I*ZA(I3-1» 

T7=T7*S:S0 

ZB(I3*ll=-T7*(T25»2. ) *Z9  ( 1 3) 

ZA(I3*1)=-T7*S3+(T25+3.)*ZA(I3> 

MF=MF*1 

GO  TO  ( 5397,5399,5401, 5413, 5405,5407,5409,5411,5413,5415,5 131) ,HF 
F 

CONTINUE 

V8(4)=-(Z6(5>-. 16444  44444  44 «44*ZB ( 6 I ♦ . 23148  14814  814  6 IE-2  * ZB( 7 
))  *.41666  66666  66667E-1 

VA(3)=(ZA(5)-.2C714  28571  42 8 57* Z A ( 6 )♦. OC 625* ZA ( 7) )*. 37 0 37  03703 
7C3  70E-1 
GO  TO  5395 
CONTINUE 

VB(5)=-(ZB(6)-. 24570  10582  0 1 355*ZB (7 ) ♦ . 1 1342  59259  25926E-1*ZQ(8) 
)-.  60281  €358:  2469 ' E-4* Z 3 (9 )> * . 3 733 3 33333  3J333E-1 
VA (4) =• 312f  3 T3rJ',  30 30 3E-1* Z A (6 > 15 . 97  00  1 7 € 3 66S4.E-2  *ZA  (7  ) 

♦ .52336  8t>067  r,194CE-3*ZA  ( 8 >-.64300  41152  263  3 7E-5*  Z A ( 5 ) 

GO  TO  5395 
CONTINUE 


f 


W8«6»=-  ( Z2(7  31303  851.67  5 2 634»  ZB  ( 8 ) ♦ . 2 38 ' 6 21693  1 21 E9E-1 *ZB( 9 

1 I-.47C67  50123  45679E- 3* 7 B ( 10 > ♦ . 1 2C 56  32716  04938E-5*ZB  ( 111) * 

2 .27777  77777  77778E-1 

VA(5)=.25C41  02564  1 0 25 bE- 1 * Z A < 7 1 - . 87 89 2 14900  32601E-2 *ZA (6 > ♦ 

1 .7892*.  16225  74956E-3*  ZA  ( <5 ) - .21568  19517  93C  63E-4*Z A ( 1 J ) 

2 ♦ 2 A < 111*.  13395  91906  72154E-6 
GO  TO  5395 

5403  CONTINUE 

V A ( 6 ) =♦ . 2222 2 22222  22Z22E-1 • 2 # 18 > - . 880 26  1 71  25  95547E- 2*ZA (9 )♦ 

1 . 10129  13208  6804CE-2*ZA <1C )-. 42538  29446  73G60E-4* 2A  (11  )♦ 

2 .62437  74137  l 3 ?2 3E-6* Zt (1 2 ) -. 22326  53177  86923E-e*ZA ( 13) 

VB(7)=-.2o8C9  5^381  95258  E-l  * Z B(  8 > ♦ . 8 81  8 5 926  72  75 934 E-2 * ZB ( 9 > - 

1 .90642  ..9258  2 32  79E -3  * ZOi  ( 1 C I ♦ . 315  1 1 18214  77  5 E 2E-4*  ZB  (1 1 ) - 

2 . 32596  73639  6 8 90  7 E-6»  ZB  (1  2 )♦. 46513  60787  22756E- «*ZB  (13) 

GO  TO  539 i 


5405  CONTINUE 

VA ( 7 ) =♦ « ,962  7 34313  7254 SE-1 * Z A ( 9 > - . 8 6 8 36  00350  Z6702E-2*ZA ( i: )♦ 
.11961  47:2  7 0 3957E-2»ZA(UI-.  66655  77805  7 1 59EE-4*  2A  ( 12  ) ♦ 
.15617  53570  8 8423 E-5* Z A (1 3 ) -. 137 99  92297  5591 7E-7*ZA (14 ) * 

• 3 1 0 O'  9 C7191  48504E-10*ZA(15) 

V8 ( 8 ) =- • 20  33  3 33333  33333E -1  • Z B(  9)  ♦ . 87544  99958  2i.  3 E6E-2*  ZB  ( 1 C>  - 
. 11192  440  92  79237E-2*ZB(1 1 ) ♦ . 5 3 3 8 41632  21  7596-  4*  ZB  (12)- 
. 1036C  16990  95848E-5*ZC(12)+. 63255  50670  62948E-8’ Z9 ( 14 ) - 
.55373  34270  5 C90G E-ll* ZB< 1 5 ) 

GO  TO  5395 
5407  CONTINUE 

V A ( 8 ) =♦  . 1754  3 85964  91228E-1 * ZA(  10 )-.  84985  14229  1 398 IE- 2*Z A ( 11)  ♦ 
.1  3448  6C333  5 0 672E -2  * Z A ( 1 2 ) - . 92 0 6 1 86407  7C  1 2 7E -4* ZA ( 1 3 ) ♦ 
.29398  2 3b0  7 35374E-5*ZA (1 4 ) -. 42707  63589  500  73t-7*ZA  ( 15  ) ♦ 
.24541  46548  6 8 95 9 E-9 * Z i ( . 6 ) - • 369 1 5 5619C  3 3 93 2E - 12  * ' A ( 1 7 ) 
VB(9>=-. If  518  51  851  851  8 5E-1  * ZC < 1 0 ) ♦ . 3 5 566  42554  975 8 9E-2* ZB ( 11) - 
.12744  2 C6  4 6 38934E-2*ZO(l2)». 79284  t>2l  52  6 8 942E-4* ZB ( 1 3 ) - 
.2198  3 32488  9550  4 E-5  * ZU  1 4 ) ♦ . 25640  60419  0 2 54  8E  - 7*  ZB  ( 1 5 ) - 

• 10G2  6 26658  5 83 16  E-9* Z B( 1 6 )♦  .5768 0 56531  7802  .E-l 3»ZB (17  ) 

GO  TO  5395 
CONTINUE 

VA(9)=*. 15873  C 15  37  3C 1 59E-1 * Z A( 11  ) -.82821  5E153  6 731 8E -2* Z A (12 ) 

♦ .14652  C 6930  61 43 9E-2» Z A (l  3 ) -. 1 175 7 92639  58914E-3*ZA ( 14) ♦ 
.47C65  76990  79954 E-5* Zi ( 15 )-. 952r 1 29452  98324t-7*ZA (16) ♦ 
.92092  73029  }935lE-9*ZA(17 1-.36402  25970  5 5366E- 11*Z A (1 8 ) ♦ 

. 38453  7 1 l 2 1 1 868CE-14*ZA(19) 

VB(  10  ) =-.16666  66666  66667E- 1 ’ZBUl  )♦  .83926  0 1741  7 3222E-2*ZB  (12  ) 
-.14062  22299  54  5 3 E-2 • ZB ( 1 3 » ♦.10 485  97453  6 1 7 35E- 3*ZB  ( 1 4 ) - 
.37786  94453  3 39  8<.  E-5  * Z6  <1  5 ) ♦ . 657  09  96457  5 7 853l  - 7*  ZB  ( 16  ) - 
.50483  50738  8 61 2 5E-9» Z 9 (1 7 » ♦ . 1 34 1 2 65412  1 e996E - 11 *Ze ( 1 8 ) - 
.53407  93084  98167E-15*ZB( 1 9) 

GO  TO  5395 
5411  CONTINUE 

VA  ( D)  = ♦ .14492  75362  31884E- 1 *ZA ( 12) - . 8C 534  4EC<>5  241  76E-2»Z  A (13  ) 

♦ .1562  5 36761  62 32 7E-2 • Z A ( 14  I -. 142 49  827E0  4 39Q lE- 3*  ZA ( 1 5 ) ♦ 
.67947  4 *.263  92697E-5* Z A (1 6 ) -. 1751*  13319  5 771 5t -6* 74 ( 1 7 > ♦ 

. 24.07  91C70  41612E-«»ZA <1  8 ) 16238  1E309  647 1 3t - 10 »Z A ( 1 9) ♦ 

. 46225  C 21  9 3 2 7 80  8E-1 3 * Z A ( 2 0 )- . 35  60  5 .?  8 7 2 3 321.UE-1  €*ZA  (21  ) 

VB( 11)=-. 15151  51515  1515:E-1*Z0(12)».31E!5  47795  4 1 82 8E- 2*Z B (1 3 ) 
-.  15164  52840  22835E-2*Z0(14M.13O14  4 5605  3 8248E-3*  ZB  ( 1 5 ) - 


1 

2 

3 

1 

2 

3 

5409 

1 

2 

3 

4 

1 

2 

3 

4 


•57147  9941  7 62549E-5»Z8(16>4. 13156  76549  2 7232E-6*  ZB  (17  >- 
.15380  56544  74373E-8»Z6<18 >♦. 82135  07414  8C055E-11 *Zu ( 19> - 
.15509  66311  87d68E-13*ZB(2C )♦ .44506  6C924  15135E-1 7* Ze (21  ! 

GO  TO  5395 
CONTINUE 

VA(11)=4. 13333  33333  33  3 33E-1  * ZA  (1  3 ) - . 7 8228  731*9  4 7945 E -2* Z A 114 ) 
♦.16412  35305  54189E-2»Zi (15 ) -.16641  22539  53536C-3*ZA < 16 ) 4 
.91369  C 96  2 7 54997E-5»ZA (17  )-.  23427  17670  5 £ 726c- ZA  118  > 4 
.50343  42785  294C 5 E-fl »Z A (1 9 ) 491 1 8 23691  34 03  1C  - 10 »Z A ( 20 > 4 
.24317  2 716  3 7648 IE-1 Z*l A ( 21  )-. 5 1 25 1 26773  740 51E-1 5* ZA (22  ) 4 
.29671  C 7269  43426E-18*ZA(2 3) 

VB(12)=-.  13838  58883  8 8 839E-1 * ZOI1 3 ) 4 , 79379  48362  751  39E-2*ZB  114 ) 
-. 16C39  76493  891 70 E-2* Zb (1 5 ) 4 . i 5459  76095  1 1 72 96 - 3*ZU ( 16 ) - 
.79380  61C11  2862L E-5»Z8(17 »4. 22501  75615  8273SE-6*ZB ( 15 ) - 
.35525  2 6641  36041 E-fl * Zb ( 19  ) 4 . 293 51  25130  221" 9E- It * 10 ( 2 0 > - 
.11351  5291  5 C7799E-12*ZB(  21)  4.15735  01C6  7 3390 3E-1 5*Z B ( 22  ) - 
.33717  128  j 6 1 7532 E-19*ZB( 2 3 ) 

GO  TO  5395 
CONTINUE 

WA(12)  = 4. 12345  679Q  1 2345 7E-1 * ZA ( 14 ) -.75962  63771  724  8 0 E-2  * ZA  f 1 

5)  4.17248  2C355  96  1C  3 E-2*Z A ( Id ) - . 1 8 9 1C  65796  7. 8 33E-3 * Z A ( 17 ) ♦ 
.11672  2392C  56-49E-4»ZA (1 8 ) -.42273  44701  2 794 5E-6* ZA < 19 ) 4 
.91376  18615  04 362E-8* Z A (2 C ) -. 11 65 3 40093  2 53 8 7* -9* Z A ( 21 ) 4 
.84215  10565  3504 8 t- 1 2»Z A ( 2 2 ) - . 31 31  8 5261C  8558 . E-l 4* ZA  (2 3 )4 
.503e4  3r725  1441 8E-17»ZA( 24)-. 22478  98637  45020E-2  *ZA(25) 

VB(  13)  = -.  12825  51282  C 5 12 8E- 1 * Z3 1 14 ) 4 . 77 0 83  287-7  6 6910 E -2* Z B (1 5 > - 
- .16747  2633  - 2 7 663 E- 2* ZO ( 1 6 ) 4 . 1 77 91  99594  95527E-3*ZJ ( 17 » - 
.10383  84179  C58C2E-4*Z9(18  ) 4.34987  60792  2 8262E-6*Z3 ( 19)  - 
.68844  72270  71  384  E -3  • Z9  (2  0 ) 4 . 774 94  50320  60 7 I5t - 10  * Z P ( 2 1 ) - 
.472  32  73762  287 1 5E - 12* Z 8 < 2 2 ) 4 . l 360  8 26691  6 8 372E- 1 -• ZB (2 3 ) - 
.14341  1846  89323E-17*Z6(24)4.23414  67226  51C 62E-21 *ZB ( 25 ) 

GO  TO  5395 

CONTINUE 

T=l./A 

C ST  = ( ( ( ( (CST5*T-CST4)*T  *CST3 ) *T4CST2) ’T-.11956  8765  4 32 C 99E-1)*1 
4.27777  77777  77778  E-l)*T4l. 

PAP1X=(CST*(T154T17I*T5)/RT2PI 
A= ASAV 

IF  ( SC . IT .2 . ) GO  TO  5441 

QANS^PAPIX 

ANS=l.-CAhS 

GO  TO  5451 

CONTINUE 

A NS  =P  AP 1 X 

QANS=1.-ANS 

CONTINUE 

RETURN 

this  is  part  bi 


CONTINUE 
N = C 

RTX=SQRT (X ) 

T 5 = E XP  ( - x ) 

RN  = T5/(RT  X*GAHPT5) 


THIS  IS  PART  B1  SUPPLEMENT 
IF  I X.GT.16.)  GO  TO  6031 

T 1=C0  *RT  X*  <C1*RTXMC2*RTX*<C3*RTX*(CL*RTXMC5«RTX*(C6*C7»RTXI  III  I 
I 

T 3 = 00  *RT  X* 101+RTX* <02»RTX* <0 3 *RTX* I 04+RTX* 105+RT X* (06*RTX*07  Hill 

> 

0EL=IT5  M1I/T3 
GO  TO  6 C 91 
CONTINUE 
T=1  ./X 

T1=ALPH0+T*< ALPH1*T*IALPH2*T • IALPH3*T*ALPH4I ) > 

T3=8ETC»T* (B£Tl*T*(8ET2'»T*(3ET3»T*BETLll I 
DEL  = I T5  /RTX1M1  ./RTPI  ♦T1/IT3*XI1 
GO  TO  6 C 91 
CONTINUE 
N = 1 

RN=£XP{ -XI 
OEL=RN 
CONTINUE 
PNC  X = DEL 

CEE= AF*  FLOAT  INI-1. 

CONTINUE 

IF  I N.EQ.I  I GO  TO  6111 
N = N*1 

CE£=CE£*1. 

RN=IX*RN»/CEE 
PNC  X = PH  C X«-RN 
GO  TO  6 1?  i 
CONTINUE 
R = X»">N 
QANS=PNCX 
ANS=1 •- CANS 
RETURN 

THIS  IS  PART  C 

CONTINUE 
A2NM1=1 . 

A 2N= 1 . 

B2NM1=X 

B2N=X+1.-A 

AN=1. 

CONTINUE 

A2Nhl=X* A2N»AN* A2NM1 
B2NM1=X*Q2N»AN*92NH1 
AMC=A2NK1/62NM1 
AN=AN*1. 

SA2N=AN- A 

A2N= A2NMi+SA2N*A2N 
B2N=82NMlfSA2N»B2N 
AN0  = A2N/(32N 

IF  I ABS(AN:-AMC  ■ .GE.IACC*ANO  1 I GO  TO  7031 
QANS=R* ANv 


1 

j 

I 


i:-io 


o o o 


AN S*1.-QANS 
RE  TURN 

THIS  IS  PART  0 

9CU  CONTINUE 
SNP1 =1 . 

SUH*SNP1 
CEE  = A 

9031  CONTINUE 

CEE=CEE-1. 

SNP1=(SNP1*CEEI/X 

SUM=SUM*SNP1 

IF  < A6S (SNPl t.GE.ACC  I GO  TO  9C31 
Q ANS= (R*  SUMI /X 

ans=i.-qans 

RETURN 

ENO 


E-!  1 


AC  « T 1 > 


I 

I, 


1 

( 

t 

f 


I 

t 


i 

i 


SUBROUTINE  GAHC03  ( A , H , 
DIMENSION  C < 30  ) 

OATA  ( C(I), 1-3, 18  ) / 


1 

♦ .5  772  1 

5 fco4  9 

0 15  3 3 

£♦0 

2 

-.65587 

8071  5 

20254 

£♦0 

3 

-.42002 

635u  3 

40952 

E-l 

4 

♦ . 16653 

36113 

822  91 

E ♦ 0 

5 

-.42197 

73455 

5 5443 

F-l 

6 

- .96219 

7152  7 

8 7o  97 

E-2 

7 

♦.72189 

4324  6 

66310 

E-2 

8 

-.11651 

67591 

8 590  7 

E-2 

« 

9 

-.21524 

1674  1 

1 4951E 

-3, 

1 

♦.12805 

02323 

8 3116 

E-3 

* 

2 

-.20134 

85478 

0 73  9 2 

E-4 

t 

3 

-.12534 

93482 

1 4267 

E-5 

• 

4 

♦ . 1 1330 

27..31 

981  70 

£-5 

9 

5 

-.20563 

38416 

9 7761 

£-6 

* 

6 

♦ .611 60 

9513  4 

4 0 1 ,2 

£-8 

« 

7 

♦ .5C020 

9/644 

4o922 

E-8 

ACC=. 5* AC 

I = 1 NT ( A ) 

AF=  A-FLOA  T (I» 

IF  ( AF.EQ.Q  . ) GO  TO  3091 
CEE  = AF- l . 

IF  ( AF.LT. 0.5  ) C£F=AF 

PHI =CEE 

K-3 

ANK -C <K)  *CEE 
SUM  - A NK 

3011  CONTI NJE 

CEE=CEE*PHI 
K=K  ♦ 1 

AN<=C (<) »CEE 
SUM=SUM*ANK 

IF  ( ARS  (ANO.LT.  <ACC»ABS  (SUM)  ) 
IF  ( K.EQ.  18  I SO  TO  3031 
GO  TO  3311 
3031  CONTINUE 
H=SUM 

IF  ( AF.GE.0.5  t H= (1.-A+SUM1/A 
Tl= 1.  ♦ SUM 

IF  < AF.LT. 0.5  ) T1  = AF  * T 1 
T1=1./T1 

IF  I I.EQ.O  » RETURN 
AJsAF-l. 

J=0 

3051  CONTINUE 
J=JU 
AJ=  AJU. 

T 1=  AJ*T 1 

IF  ( J.LT.I  ) GO  TO  3051 
3071  CONTINUE 

GO  TC  3151 
3091  CONTINUE 
J=0 


GO  TO  3031 


mm 


J 

j 

j 


] 


| 

1 


F.-12 


